for “Animal-vectored nutrient flows over resource gradients influence the nature of local and meta-ecosystem functioning”
This notebook contains the code used to run numerical analyses of a meta-ecosystem model in which consumers move across ecosystem borders in different ways with respect to resource gradients.
In section Model Numerical Analyses, we perform numerical analyses of our model of consumer movement in a two-patch meta-ecosystem, using the set of equilibrium equations derived in Mathematica (see Mathematica notebook in this same repository; Wolfram Research Inc. (2020)). In section Changes to biomass distribution in local and meta-ecosystem, we investigate how biomass distribution changes following different types of consumer movement and how this may affect top-down influences of consumer movement on local and meta-ecosystem functions. Sections Visual Deliverables and Summary Tables contain code to reproduce the figures and tables shown in the main text and Supplementary Materials.
In this introduction, we briefly define the contrasting consumer movement types we investigate with this model (Defining movement types), discuss separation of scales in time and space (Separation of Spatial and Temporal scales), describe the state variables and parameters in the model (State variables and parameters of the model), and describe the analyses setup (Analyses setup).
Here, we focus on consumers moving over gradients of inorganic nutrients availability between two local ecosystems. Consumer movement—either gradient neutral, along-, or against-gradient—connects these two local ecosystems in a meta-ecosystem (Loreau, Mouquet, and Holt 2003). We define consumer movement as gradient neutral when it happens among ecosystems that have equal nutrient availability. Along-gradient movement takes place when consumers move following the gradient direction. Conversely, against-gradient consumer movement runs counter to the direction of the gradient. Consistent with recent reviews (Gounand et al. 2018), in the following sections we also use diffusive and non-diffusive to refer to along- and against-gradient consumer movement, respectively. In our model, we create resource gradients by varying the relative values of the inorganic nutrient input rates in the two ecosystems comprising our meta-ecosystem (parameter Ii; see Table 1 in the main text or below for a list of model parameters and state variables).
While organisms can at times cross into an ecosystem that is equally suitable for them immediately after leaving their original one, this is not often the case. Rather, organisms spend significant amounts of time searching for a new suitable ecosystem, traversing the so-called “unsuitable” matrix (Weisser and Hassell 1996). Models of organismal movement in meta-ecosystems should account for this time spent outside ecosystem borders, as it may lead to loss of biomass and nutrients from the system (Gounand et al. 2018).
We model consumer movement in our meta-ecosystem model by combining a dispersers’ pool (Weisser and Hassell 1996) with the time scales separation (TSS) mathematical technique (Otto and Day 2011). Briefly, adding a dispersers’ pool to a classic meta-ecosystem model means that consumer do not instantly travel from donor to recipient ecosystem. Instead, they move from the donor to the dispersers’ pool and then from there to the recipient (see Figure 1c in the main text for a visual representation). This allows us to separate ecosystem processes that happen locally from those that happen regionally.
Some of these ecosystem processes, like primary productivity in terrestrial ecosystems, take place over long time frames (slow processes). Others, like consumer movement, happen over short time frames (fast processes). TSS explicitly accounts for these differences by allowing fast processes to reach a quasi-equilibrium while holding slow ones constant. The quasi-equilibrium is then plugged into the equation(s) for the slow processes to investigate how these processes evolve as the fast ones stay constant (Otto and Day 2011).
The combined use of the dispersers’ pool and TSS allows us to separate local and regional dynamics in our system and to investigate the effects of consumer movement at either spatial extent at the same time. See section Model Derivation in the main text for further details on our novel approach to model consumer movement in meta-ecosystems.
The table below lists the model’s state variables and parameters, and is analogous to Table 1 in the main text. For a visual representation of the model, please see Figure 1c in the main text.
| Variable | Description | Units | Range |
|---|---|---|---|
| Ni | Nutrients stocks in patch i | g | >0 |
| Pi | Producers stocks in patch i | g | >0 |
| Ci | Consumers stocks in patch i | g | >0 |
| Q | Dispersers’ Pool | g | >0 |
| Parameter | Description | Units | Range |
|---|---|---|---|
| Ii | Inorganic nutrient input rate into patch i | g*t-1 | [0,20] |
| \(l\) | Inorganic nutrient output rate | t-1 | [0,10] |
| ui | Producer uptake rate in patch i | (g*t)-1 | [0,10] |
| ai | Consumer attack rate in patch i | (g*t)-1 | [0,10] |
| \(\epsilon\)i | Consumer assimilation efficiency in patch i | unitless | [0,1] |
| hi | Biomass loss rate from Pi | t-1 | [0,10] |
| di | Biomass loss rate from Ci | t-1 | [0,10] |
| g | Movement rate from C1 to Q | t-1 | [0,10] |
| m | Movement rate from Q to C2 | t-1 | [0,10] |
| c | Biomass loss rate from Q | t-1 | [0,10] |
We derived the model’s single feasible equilibrium in Mathematica (Wolfram Research Inc. (2020); see Model Derivation in the main text for further details).
This notebook assumes the following folder structure:
.
├── Code
├── Data
└── Results
This structure allows us to take advantage of relative paths when
loading data or saving figures. Note that . stands for the
root folder. If you wish to recompile this notebook, please organize
your folder structure accordingly.
This notebook uses the following packages:
library(renv) # project dependency (packages, mostly) management
library(tidyverse) # collection of packages improving base R
library(plotrix) # additional plotting functions
library(lhs) # use latin hypercube to create toy datasets
library(ggpubr) # extra themes for ggplot2
library(ggpol) # extra geoms for ggplot2
library(patchwork) # easily create and layout multi-panel plots
library(gt) # building tables in a tidyverse framework
library(scales) # scaling functions for ggplot2
library(khroma) # color-blind palettes for ggplot2
library(xaringanExtra) # extra tabbing functionality for rmarkdown documents
library(tictoc) # timing
library(RColorBrewer) # color palettes manipulation
library(MetBrewer) # additional color palettes
library(broom) # data wrangling
# Packages that need to be installed but do not require loading
# library(distill) # notebook templates
# library(rmarkdown) # build interactive documents in R
In the interest of reproducibility, we use package renv
to keep track of this notebook’s packages.
renv::snapshot()
* Lockfile written to '~/PhD/Chapter_3/Code/renv.lock'.
We compiled this notebook on a machine using the following version of R, operative system, and necessary packages:
utils:::print.sessionInfo(sessionInfo()[-8])
R version 4.2.0 (2022-04-22)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.4
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] broom_0.8.0 MetBrewer_0.2.0 RColorBrewer_1.1-3
[4] tictoc_1.0.1 xaringanExtra_0.6.0 khroma_1.8.0
[7] scales_1.2.0 gt_0.6.0 patchwork_1.1.1
[10] ggpol_0.0.7 ggpubr_0.4.0 lhs_1.1.5
[13] plotrix_3.8-2 forcats_0.5.1 stringr_1.4.0
[16] dplyr_1.0.9 purrr_0.3.4 readr_2.1.2
[19] tidyr_1.2.0 tibble_3.1.7 ggplot2_3.3.6
[22] tidyverse_1.3.1 renv_0.15.5
Here, we perform the analyses of the meta-ecosystem model described in section Numerical Analyses of the main text. If you are interested in the code used to produce the main text’s figures and tables, this can be found in sections Visual Deliverables and Summary Tables below.
We conduct three rounds of analyses, one each for the three types of consumer movement we are interested in: gradient-neutral, along-, and against-gradient. Each round of analysis comprises iterating the model over n = 10 000 parameter sets, running stability analyses, and graphing the results. The analyses will follow these steps:
In the code chunk below, we generate two sets of randomly-drawn
parameter values to use in all scenarios. We use function
lhs() to do so. One set includes parameters values scaled
0–10 (DATA1), while the other one includes parameter values
scaled 0–1 (DATA2).
# Below we generate random values for those parameters that will vary
# during the simulations: attack rates, efficiency, death rates, and movement
# rates
DATA1 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for(i in seq(1,100,1)){
# Use lhs function in R. 100 bins, 11 iterations - one for each parameter that
# is scaled between 0 and 10.
loop_1 <- randomLHS(100, 11)
# Use rbind to stack each iteration on top of the other to create a one column
# list of 10000 samples with equal distribution in 100 bins.
# Multiplied by 10 to get range of data from 0 to 10
DATA1 <- rbind(DATA1, data.frame(i=i, Result=loop_1*10))
# print the i value so you can track progress
# print(i)
# Close the loop
}
# This is for 2 parameters that is scaled 0 to 1.
DATA2 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for(i in seq(1,100,1)){
# Use lhs function in R. 100 bins, 2 iterations.
loop_2 <- randomLHS(100, 2)
# Use rbind to stack each iteration on top of the other to create a one column
# list of 200000 samples with equal distribution in 100 bins.
# Multiplied by 10 to get range of data from 0 to 10
DATA2 <- rbind(DATA2, data.frame(i=i, Result=loop_2))
# print the i value so you can track progress
# print(i)
# Close the loop
}
# Save these two datasets for later re-use
# write.csv(DATA1, file = "../Data/DATA1.csv", row.names = FALSE)
# write.csv(DATA2, file = "../Data/DATA2.csv", row.names = FALSE)
Alternatively, we can load the two datasets provided alongside this R
notebook, in folder ../Data/. These were similarly
generated using function lhs(), scaled either 0–10 or 0–1,
and then used to produce the results and figures presented in the
paper.
Loading the provided datasets is the default approach. To generate new randomly drawn parameter values, set
eval=FALSEin the chunk below andeval=TRUEin the code chunk above.
We begin our numerical analyses of our model with the control scenario, in which consumers move between two ecosystems that have equal environmental fertility. Consumer movement in this scenario is gradient-neutral, as the environment is homogeneous in terms of resource availability. The results of this first set of model iterations produce the values used at the denominator of the response ratio (LRR) used to capture changes in local and meta-ecosystem function values (see section Visual Deliverables below). Below, we describe each step of the analyses in detail: later analyses will not describe these steps again.
First, we allocate an empty dataframe PRED to contain
the results of each iteration of the model. This dataframe will contain
the parameter values and the resulting local and meta-ecosystem function
values.
# Allocate an empty data frame to save the data in, then allocate some rows and
# columns. We will replace the current numbers with results below
# Data we are calculating are (bio)mass for soil nutrients (N1, N2), producers
# (P1, P2), consumers (C1, C2), disperses (Q)
PRED <- NULL
PRED <-rbind(PRED, data.frame(TIME=seq(0,100,1), I1 = seq(0,100,1), I2 = seq(0,100,1),
l = seq(0,100,1), u1 = seq(0,100,1), u2 = seq(0,100,1),
a1 = seq(0,100,1), a2 = seq(0,100,1), h1 = seq(0,100,1),
h2 = seq(0,100,1), d1 = seq(0,100,1), d2 = seq(0,100,1),
g = seq(0,100,1), m = seq(0,100,1), c = seq(0,100,1),
e1 = seq(0,100,1), e2 = seq(0,100,1), N1 = seq(0,100,1),
P1 = seq(0,100,1), C1 = seq(0,100,1), N2 = seq(0,100,1),
P2 = seq(0,100,1), C2 = seq(0,100,1), Q = seq(0,100,1),
FLUX_P1 = seq(0,100,1), FLUX_C1 = seq(0,100,1),
FLUX_P2 = seq(0,100,1), FLUX_C2 = seq(0,100,1),
FLUX_Eco_1 = seq(0,100,1), FLUX_Eco_2 = seq(0,100,1),
FLUX_Q = seq(0,100,1), PROD_P1 = seq(0,100,1),
PROD_C1 = seq(0,100,1), PROD_P2 = seq(0,100,1),
PROD_C2 = seq(0,100,1)))
Now we run the model. We first set the values for both nutrient input
rates in ecosystem 1 (I1 in the code chunks) and ecosystem
2 (I2) at 10. Recall that this is the control scenario with
equal environmental fertility. We also set the value of the inorganic
nutrient leaching rate (l) at 0.1.
The code chunk below runs the model using a for loop
that, over 10 000 time steps, iterates the model’s single feasible
equilibrium to calculate:
N1, N2), Producers
(P1, P2), and Consumers (C1,
C2)Q) stockFLUX_X,
where X is the compartment)PROD_X,
where X is the trophic compartment)in both local ecosystems. At each time step, we draw a random value
for each parameter from either the DATA1 or
DATA2 dataset—depending on whether the parameter is scaled
0–10 or 0–1, respectively. Then, we plug these values into the model’s
equilibria and ecosystem function formulas derived in Mathematica (Wolfram Research Inc. 2020), and finally
store both the randomly drawn parameter values and ecosystem function
values in PRED.
# Here we set values for parameters that we will not vary during the simulations
# Nutrients input rate in each patch
I1 = 10
I2 = 10
# Patch leaching rate/soil loss rate
l = 0.1
# We want to start simulations at 0
i1 = 0
for (i in seq(0,9999,1)){
# go to the next row
i1 = i1 + 1
# Let's sample the parameter values from DATA. We will sample these i times as
# defined by the for loop above
# producers uptake rates
u1 = DATA1[i1,2] # in ecosystem 1
u2 = DATA1[i1,3] # in ecosystem 2
# consumers attack rates
a1 = DATA1[i1,4] # in ecosystem 1
a2 = DATA1[i1,5] # in ecosystem 2
# death rates (= recycling rate)
h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystems 2
d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
c = DATA1[i1,10] # death rate in the disperser's pool Q
# movement rate to the Dispersers' pool Q
g = DATA1[i1,11]
# movement rate from the Dispersers' pool Q
m = DATA1[i1,12]
# assimilation efficiency of consumers
e1 = DATA2[i1,2] # in ecosystem 1
e2 = DATA2[i1,3] # in ecosystem 2
PRED[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
# Nutrient stock in ecosystem 1
((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Primary Producers biomass in ecosystem 1
(d1 + g)/(a1*e1),
# Consumers biomass in ecosystem 1
(e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Nutrient stock in ecosystem 2
(-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
# Primary Producers biomass in ecosystem 2
(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
# Consumer biomass in ecosystem 2
((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
# Consumer biomass in the Dispersers' Pool (Q)
(e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
# Nutrient Flux for Primary Producers in Ecosystem 1
((d1 + g)*h1)/(a1*e1),
# Nutrient Flux for Consumers in Ecosystem 1
(d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Nutrient Flux for Primary Producers in Ecosystem 2
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
# Nutrient Flux for Consumers in Ecosystem 2
(d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
# Nutrient Flux for Ecosystem 1
((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Nutrient Flux for Ecosystem 2
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
# Nutrient Flux in the Dispersers' Pool (Q)
(c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
# Primary Productivity in Ecosystem 1
((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
# Consumer (=Secondary) Productivity in Ecosystem 1
(e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Primary Productivity in Ecosystem 2
((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)),
# Consumer (=Secondary) Productivity in Ecosystem 2
(e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
)
}
To investigate any influence of consumer movement on nutrient cycling
at regional spatial scales, we calculate local and meta-ecosystem total
biomass (Ntot,Ptot, Ctot),
meta-ecosystem nutrient flux (MetaEcoFlux) and trophic
compartment nutrient flux (FLUX_Ptot,
FLUX_Ctot) and productivity (PROD_Ptot,
PROD_Ctot). Furthermore, we manually calculate each
ecosystem’s own total recycling flux as the sum of the flux of its 2
biotic compartments, as a check on the Mathematica-derived formulas used
above.
# calculate meta-ecosystem biomass stocks for each trophic compartment,
# a.k.a., meta-ecosystem biomass
PRED$Ntot <- PRED$N1+PRED$N2
PRED$Ptot <- PRED$P1+PRED$P2
PRED$Ctot <- PRED$C1+PRED$C2
# now calculate recycling flux for the local ecosystem manually to double check
# the formulae used above are correct
PRED$FLUX_Eco1_check <- PRED$FLUX_P1 + PRED$FLUX_C1
PRED$FLUX_Eco2_check <- PRED$FLUX_P2 + PRED$FLUX_C2
# and the meta-ecosystem recycling flux
PRED$FLUX_Ptot <- PRED$FLUX_P1 + PRED$FLUX_P2
PRED$FLUX_Ctot <- PRED$FLUX_C1 + PRED$FLUX_C2
PRED$MetaEcoFlux <- PRED$FLUX_Eco_1 + PRED$FLUX_Eco_2 - PRED$FLUX_Q
# finally, calculate the production for each compartment at meta-ecosystem level
PRED <- dplyr::mutate(PRED, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PRED <- dplyr::mutate(PRED, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")
Parameter sets that produce values of the model’s state variable that are \(\leq 0\) lack biological meaning. Accordingly, in the code chunk below, we flag these parameter sets for later exclusion from the analysis.
# Equilibrium stocks less than or equal to 0 make no biological sense, so we
# will flag them in the PRED dataset.
PRED <- PRED %>% dplyr::mutate(., biosense = ifelse(N1 > 0 &
P1 > 0 &
C1 > 0 &
N2 > 0 &
P2 > 0 &
C2 > 0 &
Q > 0,
"yes",
"no"),
.after = MetaEcoFlux) %>%
dplyr::mutate_at(vars(biosense), factor)
We perform the stability analyses on the model by fitting the same sets of parameter values and equilibrium state variable values used above to the partial derivatives extracted from the model’s Jacobian matrix evaluated in Mathematica (Wolfram Research Inc. 2020). We evaluate stability of each parameter set by checking whether the real part of the leading eigenvalue of the Jacobian is positive or negative (Otto and Day 2011). If the real part of the Jacobian’s leading eigenvalue is negative, the equilibrium is stable (Otto and Day 2011).
# create an empty dataframe
MathStab <- NULL
for (i in 1:nrow(PRED)) {
# fetch parameter values from the equilibrium simulations df
I1 = PRED$I1[i]
I2 = PRED$I2[i]
l = PRED$l[i]
u1 = PRED$u1[i]
u2 = PRED$u2[i]
a1 = PRED$a1[i]
a2 = PRED$a2[i]
h1 = PRED$h1[i]
h2 = PRED$h2[i]
d1 = PRED$d1[i]
d2 = PRED$d2[i]
g = PRED$g[i]
m = PRED$m[i]
c = PRED$c[i]
e1 = PRED$e1[i]
e2 = PRED$e2[i]
# Put in the solutions to the equilibrium values here:
dN1 = PRED$N1[i]
dP1 = PRED$P1[i]
dC1 = PRED$C1[i]
dN2 = PRED$N2[i]
dP2 = PRED$P2[i]
dC2 = PRED$C2[i]
# create the Jacobian from Mathematica
# NOTE: in transposing from Mathematica, the Jacobian is here assembled by row
Jacob <- rbind(c(-l-dP1*u1, h1-dN1*u1, d1, 0, 0, 0),
c(dP1*u1, -a1*dC1-h1+dN1*u1, -a1*dP1, 0, 0, 0),
c(0, a1*dC1*e1, -d1-g+a1*e1*dP1, 0, 0, 0),
c(0, 0, 0, -l-dP2*u2, h2-dN2*u2, d2),
c(0, 0, 0, dP2*u2, -a2*dC2-h2+dN2*u2, -a2*dP2),
c(0, 0, (g*m)/(c+m), 0, a2*dC2*e2, -d2+a2*e2*dP2))
# Check to see if the real part of the leading eigenvalue is less than 0 or not - if it is than the system is stable.
if (max(Re(base::eigen(Jacob)$values)) < 0 ){
stable <- "stable"
} else {
stable <- "unstable"
}
MathStab <- rbind(MathStab, data.frame(TIME=i, I1, I2, l, u1, u2, a1, a2, h1,
h2, d1, d2, g, m, c, e1, e2, dN1, dP1,
dC1, dN2, dP2, dC2,
EiV1 = Re(eigen(Jacob)$values[1]),
EiV2 = Re(eigen(Jacob)$values[2]),
EiV3 = Re(eigen(Jacob)$values[3]),
EiV4 = Re(eigen(Jacob)$values[4]),
EiV5 = Re(eigen(Jacob)$values[5]),
EiV6 = Re(eigen(Jacob)$values[6]),
maxEv = max(Re(base::eigen(Jacob)$values)),
stable = stable,
biosense = PRED$biosense[i]))
}
MathStab$stable <- as.factor(MathStab$stable)
# separate unstable equilibria to work with later
MathStabUS <- subset(MathStab, MathStab$stable == "unstable")
According to the stability analyses run above, about 1.55% of the parameter sets produced unstable results (i.e., 155 out of 10000 iterations).
Now, let’s check if the parameter sets that produce unstable equilibria are also those that produce equilibria lacking biological meaning—i.e., those where state variables values at equilibrium are \(\leq 0\).
biononsense <- subset(PRED, PRED$biosense == "no")
identical(as.numeric(MathStabUS[ , "TIME"]), biononsense[ , "TIME"])
[1] TRUE
It appears that they are. Hence, let’s add the stable/unstable
information to PRED.
Now, let’s exclude these unstable parameter sets from the analyses below and from our results.
We are also going to sample the dataset PRED, that
contains the results of our model’s iterations, for 1000 random
iterations results. We will store these in PRED_1k and use
them in Figures 1, 2,
3, 4, 5 below
to produce some preliminary graphs of our results.
PREDpos <- filter(PRED, PRED$biosense == "yes" & PRED$stable == "stable")
# sample PREDpos to only use 1000 random simulation results
PRED_sample1000 <- PREDpos[sample(nrow(PREDpos), size = 1000), ]
PRED_1k <- droplevels(PRED_sample1000)
Before proceeding further, let’s save the PRED object to
disk, for ease of use in the future.
write.csv(PRED, file = "../Results/PRED.csv", row.names = FALSE)
We now investigate how the presence of unidirectional,
gradient-neutral consumer movement influences the meta-ecosystem’s
behavior using box-plots produced using PRED_1k. Note that
in the figures below, “Donor” always indicates ecosystem 1 and
“Recipient” always indicates ecosystem 2.
# pivot the dataset to longer format for boxplot plotting
PRED_biomass_long <- select(PRED_1k, N1:Q, Ntot:Ctot) %>%
pivot_longer(.,
names_to = "Compartment",
values_to = "Stock",
cols = c(N1, P1, C1, N2, P2, C2, Q,
Ntot, Ptot, Ctot)) %>%
dplyr::mutate(.,
Ecosystem = if_else(Compartment == "N1" |
Compartment == "C1" |
Compartment == "P1",
"Donor",
if_else(Compartment == "N2" |
Compartment == "P2" |
Compartment=="C2",
"Recipient",
if_else(Compartment == "Q",
"Dispersers Pool",
"Meta-ecosystem"))))
PRED_biomass_long$Compartment <- as_factor(PRED_biomass_long$Compartment)
PRED_biomass_long$Ecosystem <- as_factor(PRED_biomass_long$Ecosystem)
comparts_medians <- PRED_biomass_long %>%
dplyr::group_by(Compartment, Ecosystem) %>%
dplyr::summarise(., median = median(Stock),
.groups = "keep")
ggplot(data = PRED_biomass_long,
aes(x = Compartment, y = Stock, fill = Ecosystem, col = Ecosystem)) +
geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE,
outlier.shape = 25, jitter.shape = 21) +
geom_label(data = comparts_medians,
aes(y = median, label = round(median,2)),
col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
nudge_x = -0.18, label.padding = unit(0.15, "lines"),
label.size = 0.15) +
facet_grid(.~Ecosystem, scales = "free") +
scale_color_manual(values = met.brewer("Egypt", 4)) +
scale_fill_manual(values = met.brewer("Egypt", 4)) +
theme_pubr() +
coord_cartesian(y = c(0, 100)) +
theme(legend.position = "none")
Figure 1: Biomass stock values at equilibrium for local and meta-ecosystem scale, drawn using 1k randomly-drawn data points from the model simulations. Consumer movement from ecosystem 1 to ecosystem 2 via the dispersers’ pool results in foraging pressure of consumers being transferred from producers in ecosystem 1 to those in ecosystem 2. Primary producers in ecosystem 1 are thus released from foraging pressure while those in ecosystem 2 are depressed. In turn, this leads to differences in nutrient stocks between the two ecosystems at both the local and meta-ecosystem extents. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations (i.e., 10k iterations). Values on the y-axis are limited between [0, 100] to zoom in on the variation in the results.
The graph above shows a spatial trophic cascade (Knight et al. 2005). In ecosystem 2, the
influx of consumers from ecosystem 1 increases the consumer biomass, in
turn depressing the producers’ abundance and releasing the nutrient
stock from the producers’ trophic pressure. Conversely, in ecosystem 1,
consumers’ numbers fall steadily—due to the combined action of mortality
and emigration towards ecosystem 2—thus releasing the producers from the
trophic pressure exerted by consumers. In turn, this lets producers grow
unchecked, depressing the nutrient stocks in ecosystem 1. We can see
this as well by looking at the change of log10 biomass values
for each trophic compartment with respect to the two movement
parameters: g, the rate of consumer movement from ecosystem
1 to the dispersers’ pool Q (Figure 2),
and m, the rate of consumer movement from the dispersers’
pool Q to ecosystem 2 (Figure 3).
eco_levels <- c("N1" = "Nutrient", "P1" = "Producers", "C1" = "Consumers",
"N2" = "Nutrient", "P2" = "Producers", "C2" = "Consumers",
"Q" = "Dispersers", "Ntot" = "Nutrient",
"Ptot" = "Producers", "Ctot" = "Consumers")
# pivot the dataset to longer format for boxplot plotting
PRED_long <- pivot_longer(PRED_1k,
names_to = "Compartment",
values_to = "Stock",
cols = c(N1, P1, C1, N2, P2, C2, Q,
Ntot, Ptot, Ctot)) %>%
dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" |
Compartment == "C1" |
Compartment == "P1",
"Donor",
if_else(Compartment == "N2" |
Compartment == "P2" |
Compartment=="C2", "Recipient",
if_else(Compartment == "Q",
"Dispersers",
"Meta-ecosystem"))))
PRED_long$Compartment <- as_factor(PRED_long$Compartment)
PRED_long$Ecosystem <- as_factor(PRED_long$Ecosystem)
ggplot(PRED_long, aes(g, log10(Stock), col = Ecosystem)) +
geom_point(alpha = 0.15) +
geom_smooth(method = "lm", col = "white", alpha = 1, se = T, lwd = .4) +
facet_grid(.~Compartment) +
scale_color_manual(values = met.brewer("Egypt", 4)) +
theme_pubr() +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
text = element_text(size = 11))
Figure 2: Change in biomass in the meta-ecosystems’ trophic compartments with change in the movement rate from the donor ecosystem 1 to the dispersers’ pool (g). Each facet presents results for a trophic compartment, with color identifying local and meta-ecosystem. White lines are regression lines with 95% CI (grey shaded areas). Note the log10 scale on the y-axis. Graph produced using 1000 randomly-drawn data points from the model’s analyses.
ggplot(PRED_long, aes(m, log10(Stock), col = Ecosystem)) +
geom_point(alpha = 0.15) +
geom_smooth(method = "lm", col = "white", alpha = 1, se = T, lwd = .4) +
facet_grid(.~Compartment) +
scale_color_manual(values = met.brewer("Egypt", 4)) +
theme_pubr() +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
text = element_text(size = 11))
Figure 3: Change in biomass in the meta-ecostystem’s trophic compartments with change in movement rate from the disperser pool to the recipient ecosystem 2 (m). Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log10 scale on the y-axis. All other specifications as in Figure 2.
Finally, we can see a similar effect when looking at the change in
biomass in the two local ecosystems with change in the mortality rate of
consumers while in the dispersers’ pool (parameter c in the
model; Figure 4).
ggplot(PRED_long, aes(c, log10(Stock), col = Ecosystem)) +
geom_point(alpha = 0.15) +
geom_smooth(method = "lm", col = "white", alpha = 1, se = T, lwd = .4) +
facet_grid(.~Compartment) +
scale_color_manual(values = met.brewer("Egypt", 4)) +
theme_pubr() +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
text = element_text(size = 11))
Figure 4: Change in biomass in the meta-ecostystem’s trophic compartments with change in the consumer death rate while in the dispersers’ pool (c). This could represent predation risk, but also environmental hazards and stochastic events such as crossing a dangerous linear feature of the landscape—e.g., a busy highway. Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log10 scale on the y-axis. All other specifications as in Figure 2.
In the graphs below (Figure 5), we look at nutrient flux values in each ecosystem and in the meta-ecosystem for each iteration of the model.
First, we will extract the flux columns (see above) from the
PRED_1k dataframe into a new one which we will call
FluxPRED. We will then pivot this dataframe from wide to
long format, using function pivot_longer from package
tidyr.
FluxPRED <- PRED_1k %>%
dplyr::select(TIME:c,FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
tidyr::pivot_longer(., names_to = "Compartment",
values_to = "Flux",
cols = c(FLUX_P1:MetaEcoFlux)) %>%
dplyr::mutate(., Scale = if_else(Compartment == "FLUX_P1" |
Compartment == "FLUX_C1" |
Compartment == "FLUX_Eco_1",
"Donor",
if_else(Compartment == "FLUX_P2" |
Compartment == "FLUX_C2" |
Compartment == "FLUX_Eco_2",
"Recipient",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))
Now, we produce the box-and-whiskers plot of recycling flux in local and meta-ecosystem.
FluxPRED$Scale <- as.factor(FluxPRED$Scale)
FluxPRED$Compartment <- as.factor(FluxPRED$Compartment)
comparts_medians <- FluxPRED %>%
dplyr::group_by(Compartment, Scale) %>%
dplyr::summarise(., median = median(Flux), .groups = "keep")
FluxPRED %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment,
c("FLUX_P1", "FLUX_C1", "FLUX_Eco_1",
"FLUX_P2", "FLUX_C2", "FLUX_Eco_2",
"MetaEcoFlux"))) %>%
dplyr::mutate(., Scale = fct_relevel(Scale,
c("Donor",
"Recipient",
"Meta-ecosystem"))) %>%
ggplot(., aes(x = Compartment, y = Flux)) +
geom_boxjitter(outlier.intersect = TRUE, alpha = 0.25, outlier.shape = 25, jitter.shape = 21) +
geom_label(data = comparts_medians, aes(y = median, label = round(median,2)),
col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
nudge_x = -0.18, label.padding = unit(0.15, "lines"),
label.size = 0.15) +
facet_grid(.~Scale, scales = "free") +
scale_x_discrete(labels = c("FLUX_C1" = "C1", "FLUX_P1" = "P1",
"FLUX_Ptot" = "Producers", "FLUX_C2" = "C2",
"FLUX_P2" = "P2", "FLUX_Ctot" = "Consumers",
"FLUX_Eco_1" = "Donor", "FLUX_Eco_2" = "Recipient",
"MetaEcoFlux" = "Meta-ecosystem")) +
theme_pubr() +
coord_cartesian(y = c(0, 500)) +
theme(legend.position = "none", axis.title.x = element_blank())
Figure 5: Nutrient flux at local and meta-ecosystem scales when consumers move from ecosystem 1 to 2. The two ecosystems have the same fertility—i.e., basal inorganic input values are the same. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations. The y-axis is constrained between [0, 500] to show the variation in the results. Graph produced using 1000 randomly-drawn data points from the model’s analyses.
Here, we investigate the influence of consumers moving from ecosystem
1 to ecosystem 2 when there is a gradient of nutrient availability
between the two ecosystems. We will investigate two scenarios: higher
environmental fertility in the donor ecosystem 1 (i.e., I1
>> I2) and higher environmental fertility in the
recipient ecosystem 2 (i.e., I1 << I2).
Note that the environmental leaching rate of each ecosystem—i.e., the
rate at which nutrients leave local ecosystems independently of trophic
processes—is invariant and equal between the two ecosystems
(l = 0.1).
The analysis will follow the same steps as before: draw parameter
values from our DATA1 and DATA2 data objects,
simulate model behavior, run stability analyses, remove equilibria that
are unstable and/or lack biological sense, and graphically inspect and
interpret the results.
We begin by simulating a meta-ecosystem where the donor ecosystem 1 has higher environmental fertility than the recipient ecosystem 2. Consumers move along-gradient, in a diffusive way that resembles how organismal movement has been represented in previous studies (e.g., Gravel et al. 2010).
# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 > I2
I1=18
I2=2
# Leaching does not change as is still equal across ecosystems
l = 0.1
# Allocate an empty data frame to save the data in
PREDI2 <- NULL
# Allocate some rows and columns. We will replace the current numbers with results below
# Data we are calculating are (bio)mass for soil nutrients (N1, N2), producers (P1, P2), consumers (C1,C2), disperses (Q)
PREDI2 <-rbind(PREDI2, data.frame(TIME=seq(0,100,1),
I1 = I1, I2 = I2, l = l,
u1 = seq(0,100,1), u2 = seq(0,100,1),
a1 = seq(0,100,1), a2 = seq(0,100,1),
h1 = seq(0,100,1), h2 = seq(0,100,1),
d1 = seq(0,100,1), d2 = seq(0,100,1),
g = seq(0,100,1), m = seq(0,100,1),
c = seq(0,100,1),
e1 = seq(0,100,1), e2 = seq(0,100,1),
N1=seq(0,100,1), P1=seq(0,100,1),
C1=seq(0,100,1), N2=seq(0,100,1),
P2=seq(0,100,1), C2=seq(0,100,1),
Q = seq(0,100,1),
FLUX_P1 = seq(0,100,1),
FLUX_C1 = seq(0,100,1),
FLUX_P2 = seq(0,100,1),
FLUX_C2 = seq(0,100,1),
FLUX_Eco_1 = seq(0,100,1),
FLUX_Eco_2 = seq(0,100,1),
FLUX_Q = seq(0,100,1),
PROD_P1 = seq(0,100,1),
PROD_C1 = seq(0,100,1),
PROD_P2 = seq(0,100,1),
PROD_C2 = seq(0,100,1)))
# We want to start simulations at 0
i1 = 0
for (i in seq(0,9999,1)){
# go to the next row
i1 = i1 + 1
# Let's sample the parameter values from DATA. We will sample these x times as defined by the i for loop above
# producers uptake rates
u1 = DATA1[i1,2] # in ecosystem 1
u2 = DATA1[i1,3] # in ecosystem 2
# consumers attack rates
a1 = DATA1[i1,4]
a2 = DATA1[i1,5]
# death rates (= recycling rate)
h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystem 2
d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
c = DATA1[i1,10] # death rate in the disperser's pool Q
# movement rate to the Disperser's pool Q
g = DATA1[i1,11]
# movement rate from the Disperser's pool Q
m = DATA1[i1,12]
# efficiency of consumers
e1 = DATA2[i1,2] # in ecosystem 1
e2 = DATA2[i1,3] # in ecosystem 2
PREDI2[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(d1 + g)/(a1*e1),
(e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
(e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
((d1 + g)*h1)/(a1*e1),
(d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
(d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
(c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
(e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l -I1*u1)))*u2)),
(e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
)
}
# calculate meta-ecosystem biomass stocks for each trophic compartment,
# a.k.a., meta-ecosystem productivity
PREDI2$Ntot <- PREDI2$N1 + PREDI2$N2
PREDI2$Ptot <- PREDI2$P1 + PREDI2$P2
PREDI2$Ctot <- PREDI2$C1 + PREDI2$C2
# now calculate recycling flux for the local ecosystem
PREDI2$FLUX_Eco1_check <- PREDI2$FLUX_P1 + PREDI2$FLUX_C1
PREDI2$FLUX_Eco2_check <- PREDI2$FLUX_P2 + PREDI2$FLUX_C2
# and the meta-ecosystem recycling flux
PREDI2$FLUX_Ptot <- PREDI2$FLUX_P1 + PREDI2$FLUX_P2
PREDI2$FLUX_Ctot <- PREDI2$FLUX_C1 + PREDI2$FLUX_C2
PREDI2$MetaEcoFlux <- PREDI2$FLUX_Eco_1 + PREDI2$FLUX_Eco_2 - PREDI2$FLUX_Q
# Finally, calculate the production for each compartment at the meta-ecosystem scale
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")
# Equilibrium stocks less than or equal to 0 make no biological sense, so we
# will flag them in the PREDI2 dataset.
PREDI2 <- PREDI2 %>%
dplyr::mutate(., biosense = ifelse(N1 > 0 &
P1 > 0 &
C1 > 0 &
N2 > 0 &
P2 > 0 &
C2 > 0 &
Q > 0,
"yes",
"no"),
.after = MetaEcoFlux) %>%
dplyr::mutate_at(vars(biosense), factor)
Here we perform stability analyses for the model in this scenario
where I1 >> I2. The rationale and code
are the same as above.
# create an empty dataframe
MathStabI2 <- NULL
for (i in 1:nrow(PREDI2)) {
# fetch parameter values from the equilibrium simulations df
I1 = PREDI2$I1[i]
I2 = PREDI2$I2[i]
l = PREDI2$l[i]
u1 = PREDI2$u1[i]
u2 = PREDI2$u2[i]
a1 = PREDI2$a1[i]
a2 = PREDI2$a2[i]
h1 = PREDI2$h1[i]
h2 = PREDI2$h2[i]
d1 = PREDI2$d1[i]
d2 = PREDI2$d2[i]
g = PREDI2$g[i]
m = PREDI2$m[i]
c = PREDI2$c[i]
e1 = PREDI2$e1[i]
e2 = PREDI2$e2[i]
# Put in the solutions to the equilibrium values here:
dN1 = PREDI2$N1[i]
dP1 = PREDI2$P1[i]
dC1 = PREDI2$C1[i]
dN2 = PREDI2$N2[i]
dP2 = PREDI2$P2[i]
dC2 = PREDI2$C2[i]
# create the Jacobian from Mathematica
# NOTE: in transposing from Mathematica, the Jacobian is here assembled by row
Jacob <- rbind(c(-l-dP1*u1, h1-dN1*u1, d1, 0, 0, 0),
c(dP1*u1, -a1*dC1-h1+dN1*u1, -a1*dP1, 0, 0, 0),
c(0, a1*dC1*e1, -d1-g+a1*e1*dP1, 0, 0, 0),
c(0, 0, 0, -l-dP2*u2, h2-dN2*u2, d2),
c(0, 0, 0, dP2*u2, -a2*dC2-h2+dN2*u2, -a2*dP2),
c(0, 0, (g*m)/(c+m), 0, a2*dC2*e2, -d2+a2*e2*dP2))
# Check to see if the real part of the leading eigenvalue is less than 0 or not - if it is than the system is stable.
if (max(Re(base::eigen(Jacob)$values)) < 0 ){
stable <- "stable"
} else {
stable <- "unstable"
}
MathStabI2 <- rbind(MathStabI2,
data.frame(TIME=i, I1 , I2, l, u1, u2, a1, a2, h1, h2,
d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
dP2, dC2, EiV1=Re(eigen(Jacob)$values[1]),
EiV2=Re(eigen(Jacob)$values[2]),
EiV3=Re(eigen(Jacob)$values[3]),
EiV4=Re(eigen(Jacob)$values[4]),
EiV5=Re(eigen(Jacob)$values[5]),
EiV6=Re(eigen(Jacob)$values[6]),
maxEv=max(Re(base::eigen(Jacob)$values)),
stable=stable, biosense=PREDI2$biosense[i]))
}
MathStabI2$stable <- as.factor(MathStabI2$stable)
# separate unstable equilibria to work with later
MathStabI2US <- subset(MathStabI2, MathStabI2$stable == "unstable")
3.24% of the stability analyses runs produced unstable results (i.e., 324 out of 10000 iterations). Now we check if these are the same ones that produce equilibria where state variables are < 0 at equilibrium (i.e., not biologicaly plausible).
I2biononsense <- subset(PREDI2, PREDI2$biosense == "no")
identical(as.numeric(MathStabI2US[ , "TIME"]), I2biononsense[ , "TIME"])
[1] TRUE
Since it appears that they are, let’s add the stable/unstable
information to PREDI2.
Now, let’s exclude the unstable parameter sets from the analyses
below and from our results. Then, we will sample 1000 iterations out of
PREDI2 to use in later figures.
PREDI2pos <- filter(PREDI2, PREDI2$biosense == "yes" & PREDI2$stable == "stable")
# sample PREDpos to only use 1000 random simulation results
PREDI2_sample1000 <- PREDI2pos[sample(nrow(PREDI2pos), size = 1000), ]
PREDI2_1k <- droplevels(PREDI2_sample1000)
# pivot the dataframe
PREDI2_biomass_long <- pivot_longer(PREDI2_1k,
names_to = "Compartment",
values_to = "Stock",
cols = c(N1, P1, C1, N2, P2, C2, Q,
Ntot, Ptot, Ctot)) %>%
dplyr::mutate(.,
Ecosystem = if_else(Compartment == "N1" |
Compartment == "P1" |
Compartment == "C1", "Donor",
if_else(Compartment == "N2" |
Compartment == "P2" |
Compartment=="C2",
"Recipient",
if_else(Compartment =="Q",
"Dispersers' Pool",
"Meta-ecosystem")))) %>%
dplyr::mutate(.,
Ecosystem = fct_relevel(Ecosystem,
"Donor",
"Recipient",
"Meta-ecosystem")) %>%
dplyr::filter(.,
Ecosystem != "Dispersers' Pool")
PREDI2_biomass_long$Compartment <- as_factor(PREDI2_biomass_long$Compartment)
PREDI2_biomass_long$Ecosystem <- as_factor(PREDI2_biomass_long$Ecosystem)
comparts_medians <- PREDI2_biomass_long %>%
dplyr::group_by(Compartment, Ecosystem) %>%
dplyr::summarise(., median = median(Stock), .groups = "keep")
I2bm <- ggplot(data=PREDI2_biomass_long,
aes(x = Compartment, y = Stock, fill = Ecosystem, col = Ecosystem)) +
geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE, outlier.shape = 25, jitter.shape = 21) +
geom_label(data = comparts_medians,
aes(y = median, label = round(median,2)),
col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
nudge_x = -0.18, label.padding = unit(0.15, "lines"),
label.size = 0.15) +
facet_grid(.~Ecosystem, scales = "free") +
scale_color_manual(values = met.brewer("Egypt", 4)) +
scale_fill_manual(values = met.brewer("Egypt", 4)) +
ggtitle("Organic Biomass and Nutrients Stock") +
theme_pubr() +
coord_cartesian(y = c(0, 100)) +
theme(legend.position = "none")
I2bm
Figure 6: Effects of along-gradient, diffusive consumers movement in a meta-ecosystem. Compared to the control scenario, the release of primary producers in ecosystem 1 is larger following movement of consumers out of this ecosystem (confront this figure with Figure 1). Environmental fertility values: I1 = 18, I2 = 2. All other specifications as in Figure 1.
FluxPREDI2 <- PREDI2_1k %>%
dplyr::select(.,
TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
tidyr::pivot_longer(.,
names_to = "Compartment",
values_to = "Flux",
cols = c(FLUX_P1:MetaEcoFlux)) %>%
dplyr::mutate(.,
Scale = if_else(Compartment == "FLUX_P1"
| Compartment == "FLUX_C1"
| Compartment == "FLUX_Eco_1",
"Donor", if_else(Compartment == "FLUX_P2" |
Compartment == "FLUX_C2" |
Compartment == "FLUX_Eco_2",
"Recipient",
"Meta-ecosystem"))) %>%
dplyr::mutate(.,
Scale = fct_relevel(Scale,
"Donor",
"Recipient",
"Meta-ecosystem"))
FluxPREDI2$Scale <- as.factor(FluxPREDI2$Scale)
FluxPREDI2$Compartment <- as.factor(FluxPREDI2$Compartment)
comparts_medians <- FluxPREDI2 %>%
dplyr::group_by(Compartment, Scale) %>%
dplyr::summarise(., median = median(Flux), .groups = "keep")
I2nf <- FluxPREDI2 %>%
dplyr::mutate(.,
Compartment = fct_relevel(Compartment,
c("FLUX_P1", "FLUX_C1", "FLUX_Eco_1",
"FLUX_P2", "FLUX_C2", "FLUX_Eco_2",
"MetaEcoFlux"))) %>%
ggplot(aes(x = Compartment, y = Flux)) +
geom_boxjitter(outlier.intersect = TRUE, alpha = 0.1, outlier.shape = 25, jitter.shape = 21) +
geom_label(data = comparts_medians,
aes(y = median, label = round(median,2)),
col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
nudge_x = -0.18, label.padding = unit(0.15, "lines"),
label.size = 0.15) +
facet_grid(.~Scale, scales = "free") +
coord_cartesian(y = c(0, 300)) +
scale_x_discrete(labels = c("FLUX_C1" = "C1",
"FLUX_P1" = "P1",
"FLUX_Eco_1" = "Donor",
"FLUX_C2" = "C2",
"FLUX_P2" = "P2",
"FLUX_Eco_2" = "Recipient",
"MetaEcoFlux" = "Meta-ecosystem")) +
ggtitle("Nutrient Flux") +
theme_pubr() +
theme(legend.position = "none", axis.title.x = element_blank())
I2nf
Figure 7: With a higher fertility in the ecosystem 1 (i.e., I1 >> I2) nutrient flux in this ecosystem increases, most likely as a consequence of the release primary producers. Conversely, ecosystem 2 shows lower and less variable levels of nutrient flux. At the meta-ecosystem scale, nutrient flux is higher than in either local ecosystems but lower than in the case of equal environmental fertility across ecosystems (Figure 5). Environmental fertility values: I1 = 18, I2 = 2. All other specifications as in Figure 5.
# I2all <- I2bm + I2nf + plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") + plot_layout(ncol = 1, nrow = 2)
#ggsave(I2all, filename = "../Results/LowI2BmNf.pdf", device = "pdf", dpi = 300, width = 8, height = 6)
The code chunk below analyses the experimental scenario of higher environmental fertility in ecosystem 2. In this case, consumer movement happens counter to the resource availability gradient (i.e., non-diffusive movement; Gounand et al. (2018)).
# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 < I2
I1=2
I2=18
# leaching rate
l = 0.1
# Allocate an empty data frame to save the data in
PREDI1 <- NULL
# Allocate some rows and columns. We will replace the current numbers with results below
# Data we are calculating are (bio)mass for soil nutrients (N1, N2), producers (P1, P2), consumers (C1,C2), disperses (Q)
PREDI1 <-rbind(PREDI1,
data.frame(TIME=seq(0,100,1), I1 = I1, I2 = I2, l = l,
u1 = seq(0,100,1), u2 = seq(0,100,1), a1 = seq(0,100,1),
a2 = seq(0,100,1), h1 = seq(0,100,1), h2 = seq(0,100,1),
d1 = seq(0,100,1), d2 = seq(0,100,1), g = seq(0,100,1),
m = seq(0,100,1), c = seq(0,100,1), e1 = seq(0,100,1),
e2 = seq(0,100,1), N1=seq(0,100,1), P1=seq(0,100,1),
C1=seq(0,100,1), N2=seq(0,100,1), P2=seq(0,100,1),
C2=seq(0,100,1), Q = seq(0,100,1), FLUX_P1 = seq(0,100,1),
FLUX_C1 = seq(0,100,1), FLUX_P2 = seq(0,100,1),
FLUX_C2 = seq(0,100,1), FLUX_Eco_1 = seq(0,100,1),
FLUX_Eco_2 = seq(0,100,1), FLUX_Q = seq(0,100,1),
PROD_P1 = seq(0,100,1), PROD_C1 = seq(0,100,1),
PROD_P2 = seq(0,100,1), PROD_C2 = seq(0,100,1)))
# We want to start simulations at 0
i1 = 0
for (i in seq(0,9999,1)){
# go to the next row
i1 = i1 + 1
# Let's sample the parameter values from DATA. We will sample these x times as defined by the i for loop above
# producers uptake rates
u1 = DATA1[i1,2] # in ecosystem 1
u2 = DATA1[i1,3] # in ecosystem 2
# consumers attack rates
a1 = DATA1[i1,4]
a2 = DATA1[i1,5]
# death rates (= recycling rate)
h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystem 2
d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
c = DATA1[i1,10] # death rate in the disperser's pool Q
# movement rate to the Disperser's pool Q
g = DATA1[i1,11]
# movement rate from the Disperser's pool Q
m = DATA1[i1,12]
# efficiency of consumers
e1 = DATA2[i1,2] # in ecosystem 1
e2 = DATA2[i1,3] # in ecosystem 2
PREDI1[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(d1 + g)/(a1*e1),
(e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
(e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
((d1 + g)*h1)/(a1*e1),
(d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
(d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
(c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
(e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l -I1*u1)))*u2)),
(e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
)
}
# calculate meta-ecosystem biomass stocks for each trophic compartment,
# a.k.a., meta-ecosystem productivity
PREDI1$Ntot <- PREDI1$N1 + PREDI1$N2
PREDI1$Ptot <- PREDI1$P1 + PREDI1$P2
PREDI1$Ctot <- PREDI1$C1 + PREDI1$C2
# now calculate recycling flux for the local ecosystem
PREDI1$FLUX_Eco1_check <- PREDI1$FLUX_P1 + PREDI1$FLUX_C1
PREDI1$FLUX_Eco2_check <- PREDI1$FLUX_P2 + PREDI1$FLUX_C2
# and the meta-ecosystem recycling flux
PREDI1$FLUX_Ptot <- PREDI1$FLUX_P1 + PREDI1$FLUX_P2
PREDI1$FLUX_Ctot <- PREDI1$FLUX_C1 + PREDI1$FLUX_C2
PREDI1$MetaEcoFlux <- PREDI1$FLUX_Eco_1 + PREDI1$FLUX_Eco_2 - PREDI1$FLUX_Q
# Finally, calculate the compartment production at meta-ecosystem scale
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")
# Equilibrium stocks less than or equal to 0 make no biological sense, so we
# will remove them from the PRED dataset.
PREDI1 <- PREDI1 %>%
dplyr::mutate(., biosense = ifelse(N1 > 0 &
P1 > 0 &
C1 > 0 &
N2 > 0 &
P2 > 0 &
C2 > 0 &
Q > 0,
"yes",
"no"),
.after = MetaEcoFlux) %>%
dplyr::mutate_at(vars(biosense), factor)
Here we perform stability analyses for the scenario of higher
environmental fertility in ecosystem 2, the recipient (i.e.,
I1 << I2). The analyses follow the
established workflow (see above).
# create an empty dataframe
MathStabI1 <- NULL
for (i in 1:nrow(PREDI1)) {
# fetch parameter values from the equilibrium simulations df
I1 = PREDI1$I1[i]
I2 = PREDI1$I2[i]
l = PREDI1$l[i]
u1 = PREDI1$u1[i]
u2 = PREDI1$u2[i]
a1 = PREDI1$a1[i]
a2 = PREDI1$a2[i]
h1 = PREDI1$h1[i]
h2 = PREDI1$h2[i]
d1 = PREDI1$d1[i]
d2 = PREDI1$d2[i]
g = PREDI1$g[i]
m = PREDI1$m[i]
c = PREDI1$c[i]
e1 = PREDI1$e1[i]
e2 = PREDI1$e2[i]
# Put in the solutions to the equilibrium values here:
dN1 = PREDI1$N1[i]
dP1 = PREDI1$P1[i]
dC1 = PREDI1$C1[i]
dN2 = PREDI1$N2[i]
dP2 = PREDI1$P2[i]
dC2 = PREDI1$C2[i]
# create the Jacobian from Mathematica
# NOTE: in transposing from Mathematica, the Jacobian is here assembled by row
Jacob <- rbind(c(-l-dP1*u1, h1-dN1*u1, d1, 0, 0, 0),
c(dP1*u1, -a1*dC1-h1+dN1*u1, -a1*dP1, 0, 0, 0),
c(0, a1*dC1*e1, -d1-g+a1*e1*dP1, 0, 0, 0),
c(0, 0, 0, -l-dP2*u2, h2-dN2*u2, d2),
c(0, 0, 0, dP2*u2, -a2*dC2-h2+dN2*u2, -a2*dP2),
c(0, 0, (g*m)/(c+m), 0, a2*dC2*e2, -d2+a2*e2*dP2))
# Check to see if the real part of the leading eigenvalue is less than 0 or not - if it is than the system is stable.
if (max(Re(base::eigen(Jacob)$values)) < 0 ){
stable <- "stable"
} else {
stable <- "unstable"
}
MathStabI1 <- rbind(MathStabI1,
data.frame(TIME=i, I1 , I2, l, u1, u2, a1, a2, h1, h2,
d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
dP2, dC2, EiV1=Re(eigen(Jacob)$values[1]),
EiV2=Re(eigen(Jacob)$values[2]),
EiV3=Re(eigen(Jacob)$values[3]),
EiV4=Re(eigen(Jacob)$values[4]),
EiV5=Re(eigen(Jacob)$values[5]),
EiV6=Re(eigen(Jacob)$values[6]),
maxEv=max(Re(base::eigen(Jacob)$values)),
stable=stable, biosense=PREDI1$biosense[i]))
}
MathStabI1$stable <- as.factor(MathStabI1$stable)
# separate unstable equilibria to work with later
MathStabI1US <- subset(MathStabI1, MathStabI1$stable == "unstable")
2.85% of the stability analyses runs produced unstable results (i.e., 285 out of 10000 iterations). As always, let’s check whether these are also the parameter sets that produce state variable values at equilibrium that do not have biological meaning (i.e., <0).
I1biononsense <- subset(PREDI1, PREDI1$biosense == "no")
identical(as.numeric(MathStabI1US[ , "TIME"]), I1biononsense[ , "TIME"])
[1] TRUE
Now, since that proved to be the case, let’s add the stable/unstable
information to PREDI1.
Now, let’s exclude the unstable, biological nonsense parameter sets
from the analyses below and from our results. Then, we sample 1000
random iterations from PREDI1 and store them for later use
in figures and comparisons.
PREDI1pos <- filter(PREDI1, PREDI1$biosense == "yes" & PREDI1$stable == "stable")
# sample PREDpos to only use 1000 random simulation results
PREDI1_sample1000 <- PREDI1pos[sample(nrow(PREDI1pos), size = 1000), ]
PREDI1_1k <- droplevels(PREDI1_sample1000)
# pivot the
PREDI1_biomass_long <- pivot_longer(PREDI1_1k, names_to = "Compartment", values_to = "Stock", cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot, Ctot)) %>% dplyr::mutate(., Ecosystem = if_else(Compartment == "N1"|Compartment == "C1"|Compartment == "P1", "Donor", if_else(Compartment == "N2"|Compartment == "P2"|Compartment=="C2", "Recipient", if_else(Compartment == "Q", "Dispersers' Pool", "Meta-ecosystem")))) %>% dplyr::mutate(., Ecosystem = fct_relevel(Ecosystem, "Donor", "Recipient", "Meta-ecosystem")) %>% dplyr::filter(., Ecosystem != "Dispersers' Pool")
PREDI1_biomass_long$Compartment <- as_factor(PREDI1_biomass_long$Compartment)
PREDI1_biomass_long$Ecosystem <- as_factor(PREDI1_biomass_long$Ecosystem)
comparts_medians <- PREDI1_biomass_long %>% dplyr::group_by(Compartment, Ecosystem) %>% dplyr::summarise(., median = median(Stock), .groups = "keep")
I1bm <- ggplot(data=PREDI1_biomass_long, aes(x = Compartment, y = Stock, fill = Ecosystem, col = Ecosystem)) +
geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE, outlier.shape = 25, jitter.shape = 21) +
geom_label(data = comparts_medians, aes(y = median, label = round(median,2)), col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75, nudge_x = -0.18, label.padding = unit(0.15, "lines"), label.size = 0.15) +
facet_grid(.~Ecosystem, scales = "free") +
scale_color_manual(values = met.brewer("Egypt", 4)) +
scale_fill_manual(values = met.brewer("Egypt", 4)) +
ggtitle("Organic Biomass and Nutrient Stock") +
theme_pubr() +
coord_cartesian(y = c(0, 100)) +
theme(legend.position = "none")
I1bm
Figure 8: Movement of consumer against the environmental fertility gradient leads to a strong reduction of the nutrients stock in the donor ecosystem 1, stemming from the release of local primary producers from consumer control. In the recipient ecosystem 2, the nutrients stock increases beyond the level of enrichment seen in Figure 1, as local and immigrant consumer strongly suppress the local primary producers. All other specifications as in Figure 1. Environmental fertility values: I1 = 2, I2 = 18.
Following from the stronger trophic cascades, nutrient flux is also affected at both local and meta-ecosystem scales when consumer movement happens against the fertility gradient. Here, median nutrient flux in the recipient, and more fertile, ecosystem 2 is more than double the median flux in the case of diffusive consumer movement (Figure 7). Furthermore, it is also higher than in the case of equal local environmental fertility—i.e., lack of a fertility gradient (Figure 5). Likewise, median nutrient flux at the meta-ecosystem scale is higher than in either the along-gradient (Figure 7) and no-gradient (Figure 5) cases.
FluxPREDI1 <- PREDI1_1k %>% dplyr::select(TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>% tidyr::pivot_longer(names_to = "Compartment", values_to = "Flux", cols = c(FLUX_P1:MetaEcoFlux)) %>% dplyr::mutate(Scale = if_else(Compartment == "FLUX_P1" | Compartment == "FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment == "FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2", "Recipient", "Meta-ecosystem"))) %>% dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))
FluxPREDI1$Scale <- as.factor(FluxPREDI1$Scale)
FluxPREDI1$Compartment <- as.factor(FluxPREDI1$Compartment)
comparts_medians <- FluxPREDI1 %>% dplyr::group_by(Compartment, Scale) %>% dplyr::summarise(., median = median(Flux), .groups = "keep")
I1nf <- FluxPREDI1 %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1", "FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
ggplot(., aes(x = Compartment, y = Flux)) +
geom_boxjitter(outlier.intersect = TRUE, alpha = 0.1, outlier.shape = 25, jitter.shape = 21) +
geom_label(data = comparts_medians, aes(y = median, label = round(median,2)), col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75, nudge_x = -0.18, label.padding = unit(0.15, "lines"), label.size = 0.15) +
facet_grid(.~Scale, scales = "free") +
theme_pubr() +
coord_cartesian(y = c(0, 400)) +
scale_x_discrete(labels = c("FLUX_C1" = "C1", "FLUX_P1" = "P1", "FLUX_Eco_1" = "Donor", "FLUX_C2" = "C2", "FLUX_P2" = "P2", "FLUX_Eco_2" = "Recipient", "MetaEcoFlux" = "Meta-ecosystem")) +
ggtitle("Nutrient Flux") +
theme_pubr() +
theme(legend.position = "none", axis.title.x = element_blank())
I1nf
Figure 9: Nutrient flux more than doubles in ecosystem 2, when this is the more fertile of the two ecosystem and is also the recipient of consumer movement. In turn, this leads to an overall increase in the nutrient flux at the meta-ecosystem scale, compared to scenarios where environmental fertility is equal among local ecosystem (Figure 5) or is higher in the donor ecosystem (Figure 19). Environmental fertility values: I1 = 2, I2 = 18.
# I1all <- I1bm + I1nf + plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") + plot_layout(ncol = 1, nrow = 2)
#ggsave(I1all, filename = "../Results/LowI1BmNf.pdf", device = "pdf", dpi = 300, width = 8, height = 6)
Here, we summarize the results of the model’s stability analyses, and collect this information to produce Table S1 shown in the Supplementary Materials. As the table shows, in all three scenario tested, only a small portion of the parameter sets used to simulate the model’s behavior at equilibrium produces unstable results. As detailed above, we exclude these unstable parameter sets from further investigations.
The table is saved in folder ../Results/ as a
.tex file.
# summarise the number of stable and unstable parameter sets for each
# stability analysis run with the numerical approximation method and
# calculate the percentage of unstable parameter sets over the total
# iterations of the simulations
MStabSum <- MathStab %>%
dplyr::group_by(., stable) %>%
dplyr::summarise(., n = n(), .groups = "keep") %>%
add_column(., Fertility = "Equal fertility") %>%
pivot_wider(id_cols = Fertility,
names_from = stable,
values_from = n) %>%
dplyr::mutate(., Percent = ((unstable/nrow(MathStab))*100))
I2MStabSum <- MathStabI2 %>%
dplyr::group_by(., stable) %>%
dplyr::summarise(., n = n(), .groups = "keep") %>%
add_column(., Fertility = "Higher fertility in ecosystem 1") %>%
pivot_wider(id_cols = Fertility,
names_from = stable,
values_from = n) %>%
dplyr::mutate(., Percent = ((unstable/nrow(MathStabI2))*100))
I1MStabSum <- MathStabI1 %>%
dplyr::group_by(., stable) %>%
dplyr::summarise(., n = n(), .groups = "keep") %>%
add_column(., Fertility = "Higher fertility in ecosystem 2") %>%
pivot_wider(id_cols = Fertility,
names_from = stable,
values_from = n) %>%
dplyr::mutate(., Percent = ((unstable/nrow(MathStabI1))*100))
# Bind the objects generated above into a single dataframe and then produce a
# table that shows the number of stable and unstable parameter sets for each
# stability analysis method
SAsummary <- bind_rows(MStabSum, I2MStabSum, I1MStabSum) %>%
ungroup() %>%
gt() %>%
fmt_number(columns = 4, decimals = 2) %>%
cols_label(., Fertility = "Fertility",
stable = md("Stable"),
unstable = md("Unstable"),
Percent = md("\\%")) %>%
tab_header(title = md("**Summary of stability analyses**, showing the number of
Stable and Unstable equilibria, and the percentage of Unstable
equilibria over the total number of iterations (n = 10000).
Rows are grouped by environmental fertility gradient scenario.")) %>%
tab_style(style = list(cell_text(weight = "bold")),
locations = cells_column_labels()) %>%
tab_style(style = cell_text(align = "left", size = 18),
locations = cells_title())
# save it
# gtsave(SAsummary, filename = "../Results/StabilityAnalyses_SummaryTable.tex")
# view it
SAsummary
| Summary of stability analyses, showing the number of Stable and Unstable equilibria, and the percentage of Unstable equilibria over the total number of iterations (n = 10000). Rows are grouped by environmental fertility gradient scenario. | |||
|---|---|---|---|
| Fertility | Stable | Unstable | % |
| Equal fertility | 9845 | 155 | 1.55 |
| Higher fertility in ecosystem 1 | 9676 | 324 | 3.24 |
| Higher fertility in ecosystem 2 | 9715 | 285 | 2.85 |
Here, we investigate whether different types of consumer movement between the two local ecosystems lead to changes in the distribution of biomass in the meta-ecosystem. Biomass distribution in local and meta-ecosystems can vary from bottom-heavy to top-heavy (McCauley et al. 2018). Bottom-heavy ecosystems present the classic biomass pyramid, with a large base of inorganic nutrients and primary producers supporting an increasingly smaller biomass of primary and secondary consumers. This bottom-heavy distribution is typical of terrestrial ecosystems—albeit exceptions exist (Hatton et al. 2015; McCauley et al. 2018). Conversely, in top-heavy ecosystems, a small base of inorganic nutrients and primary producers supports a large biomass of primary and secondary consumers. This inverted biomass pyramid, long seen as an exception, is increasingly recognized as widespread in certain systems—e.g., marine ecosystems (Sandin and Zgliczynski 2015; McCauley et al. 2018; Woodson, Schramski, and Joye 2018).
To investigate biomass distribution in our model, we use the Consumer to Resource biomass ratio, calculated as:
\(\displaystyle C{:}R = \frac{Biomass_{C^{\ast}}}{Biomass_{R^{\ast}}}\)
where the asterisk \(^{\ast}\) indicates the equilibrium biomass estimates. Values of \(C{:}R < 1\) identify ecosystems with a bottom-heavy biomass pyramid, whereas when \(C{:}R > 1\) the biomass pyramid is inverted and top-heavy (McCauley et al. 2018).
We begin by calculating the \(C{:}R\) ratio for each scenario
investigated above. First, we create a new dataframe
bmDistr, where we will store the required data: fertility
conditions, movement rates values, efficiency rates values, state
variable values at equilibrium, and ecosystem function values. We will
also filter out those parameter sets that have no biological sense and
are unstable as per the stability analyses above.
Now, we calculate \(C{:}R\). For our purposes, biomass values of primary producers in the local and meta-ecosystem will be the denominator of the ratio. That is, \(R = P^*_i\).
bmDistr <- dplyr::mutate(bmDistr,
Eco1 = C1/P1,
Eco2 = C2/P2,
MetaEco = Ctot/Ptot
)
Let’s look at the density distribution of these data. First, we exclude parameter sets that produce extreme values of \(C{:}R\)—that is, \(C{:}R > 10\).
ecoCR <- filter(bmDistr, Eco1 <= 10 & Eco2 <= 10 & MetaEco <= 10)
# summary(ecoCR)
Now, we pivot the resulting dataset and plot the density distribution of \(C{:}R\) values.
fert.labs <- c("Against-gradient", "Along-gradient", "Gradient-neutral")
names(fert.labs) <- c("Against", "Along", "Equal")
scale.labs <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")
names(scale.labs) <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")
ecoCR_long <- pivot_longer(ecoCR, cols = c(Eco1, Eco2, MetaEco),
values_to = "CRratio", names_to = "Scale") %>%
mutate_at(., vars(Scale), factor) %>%
mutate(., Scale = ifelse(Scale == "Eco1",
"Ecosystem 1", ifelse(Scale == "Eco2",
"Ecosystem 2",
"Meta-ecosystem")))
ggplot(ecoCR_long, aes(x = CRratio, fill = Fertility, col = Fertility)) +
geom_density(aes(y = ..density..), alpha = 0.1) +
scale_color_highcontrast(reverse = T, name = "Consumer movement", labels = c("Against-gradient", "Along-gradient", "Gradient-neutral")) +
scale_fill_highcontrast(reverse = T, name = "Consumer movement", labels = c("Against-gradient", "Along-gradient", "Gradient-neutral")) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
theme_pubr(legend = "right", base_size = 10) +
xlab("C:R ratio")
Figure 10: Density distribution when dataset is limited to C:R < 10. Note the different scales on the y-axis.
ggsave(filename = "../Results/Density_CRratio_lessThan10.pdf", device = "pdf", dpi = 300, width = 10, height = 7)
Here we look at how the \(C{:}R\)
values change with increasing movement from ecosystem 1, the donor,
towards the dispersers’ pool Q (i.e., parameter
g in the model). The graph below uses the full dataset
created above, ecoCR_long, and fits linear models to each
combination of Fertility scenario (gradient-neutral,
along-, and against-gradient) and Scale (ecosystem 1,
ecosystem 2, and meta-ecosystem.
ecoCRg_lm <- ecoCR_long %>% nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ g, data = data)))
ecoCRg_lmOutput <- ecoCRg_lm %>% summarise(glance(model))
ecoCRg_lmOutput$g_slope <- NA
ecoCRg_lmOutput$g_se <- NA
for (i in 1:length(ecoCRg_lm[["model"]])) {
ecoCRg_lmOutput$g_slope[i] <- ecoCRg_lm$model[[i]]$coefficient[2]
ecoCRg_lmOutput$g_se[i] <- summary(ecoCRg_lm$model[[i]])$coefficient["g", "Std. Error"]
}
ecoCRg_lmOutput %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
g_slope = "g Coeff.",
g_se = "g SE") %>%
gtsave(., filename = "../Results/CR_g_lm_res.rtf")
ecoCR_long %>% ggplot(., aes(x = g, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values [0,10].") +
theme_pubr(legend = "right", base_size = 10)
Figure 11: Distribution of C:R values lower than 10 as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Darker shades of blue correspond to a higher count of C:R value in a given bin. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data. C:R values appear to follow a decreasing trend as g increases in Ecosystem 1 and in the Meta-ecosystem, whereas they appear stable or slightly increasing in Ecosystem 2.
ggsave(filename = "../Results/CRratioVg_binned.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
Now, let’s explore the relationships in a bit more detail. We do this by splitting the dataframe in two portions: \(C{:}R > 1\) and \(C{:}R < 1\). Here, we take a look at the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).
ecoCRg_lm_above1 <- ecoCR_long %>%
filter(., CRratio > 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ g, data = data)))
ecoCRg_lmOutput_above1 <- ecoCRg_lm_above1 %>% summarise(glance(model))
ecoCRg_lmOutput_above1$g_slope <- NA
ecoCRg_lmOutput_above1$g_se <- NA
for (i in 1:length(ecoCRg_lm_above1[["model"]])) {
ecoCRg_lmOutput_above1$g_slope[i] <- ecoCRg_lm_above1$model[[i]]$coefficient[2]
ecoCRg_lmOutput_above1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g", "Std. Error"]
}
ecoCRg_lmOutput_above1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
g_slope = "g Coeff.",
g_se = "g SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers'
Pool (g)"),
subtitle = md("C:R values between (1,10].")) %>%
gtsave(., filename = "../Results/CR_g_lm_above1.rtf")
ecoCR_long %>%
filter(., CRratio > 1) %>%
ggplot(., aes(x = m, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values between (1,10].") +
theme_pubr(legend = "right", base_size = 10)
Figure 12: Distribution of C:R values comprised between (1, 10] as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. As g increases, we do not observe any trend in the values of C:R. Note that C:R > 1 are generally indicative of top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVg_binned_above1.pdf", device = "pdf",
dpi = 300, width = 12, height = 12)
And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).
ecoCRg_lm_below1 <- ecoCR_long %>%
filter(., CRratio <= 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ g, data = data)))
ecoCRg_lmOutput_below1 <- ecoCRg_lm_below1 %>% summarise(glance(model))
ecoCRg_lmOutput_below1$g_slope <- NA
ecoCRg_lmOutput_below1$g_se <- NA
for (i in 1:length(ecoCRg_lm_below1[["model"]])) {
ecoCRg_lmOutput_below1$g_slope[i] <- ecoCRg_lm_below1$model[[i]]$coefficient[2]
ecoCRg_lmOutput_below1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g", "Std. Error"]
}
ecoCRg_lmOutput_below1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(., columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
g_slope = "g Coeff.",
g_se = "g SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers'
Pool (g)"),
subtitle = md("C:R values between [0,1].")) %>%
gtsave(., filename = "../Results/CR_g_lm_below1.rtf")
ecoCR_long %>%
filter(., CRratio <= 1) %>%
ggplot(., aes(x = m, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values between [0,1].") +
scale_color_manual(name = "Legend", values = c("cyan", "magenta")) +
theme_pubr(legend = "right", base_size = 10)
Figure 13: Distribution of C:R values [0,1] when movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Again, no trend is apparent save for a seeming increase in C:R values in ecosystem 2 when consumers move along-gradient (central panel). Note, however, that this trend does not scale up to the meta-ecosystem. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVg_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
Here, we investigate how the \(C{:}R\) biomass ratio changes with
increasing values of m, the movement rate of consumers
leaving the dispersers’ pool Q to enter ecosystem 2. Here
is the general relationship.
ecoCRm_lm <- ecoCR_long %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ m, data = data)))
ecoCRm_lmOutput <- ecoCRm_lm %>% summarise(glance(model))
ecoCRm_lmOutput$m_slope <- NA
ecoCRm_lmOutput$m_se <- NA
for (i in 1:length(ecoCRm_lm[["model"]])) {
ecoCRm_lmOutput$m_slope[i] <- ecoCRm_lm$model[[i]]$coefficient[2]
ecoCRm_lmOutput$m_se[i] <- summary(ecoCRm_lm$model[[i]])$coefficient["m", "Std. Error"]
}
ecoCRm_lmOutput %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,logLik:deviance, m_slope:m_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
m_slope = "m Coeff.",
m_se = "m SE") %>%
gtsave(., filename = "../Results/CR_m_lm_res.rtf")
ecoCR_long %>%
ggplot(., aes(x = m, y = CRratio, )) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values limited to [0,10].") +
scale_color_manual(name = "Legend", values = c("cyan", "magenta")) +
theme_pubr(legend = "right", base_size = 10)
Figure 14: Density distribution of C:R values [0,10] as consumer movement rate from the dispersers’ pool to ecosystem 2 (m) increases. No trend is clearly visible, aside from an apparent increase in C:R in ecosystem 2 when consumers move along-gradients (central panel). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVm_binned.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).
ecoCRm_lm_above1 <- ecoCR_long %>%
filter(., CRratio > 1) %>%
nest_by(Fertility,
Scale) %>%
mutate(model = list(lm(CRratio ~ m,
data = data)))
ecoCRm_lmOutput_above1 <- ecoCRm_lm_above1 %>% summarise(glance(model))
ecoCRm_lmOutput_above1$m_slope <- NA
ecoCRm_lmOutput_above1$m_se <- NA
for (i in 1:length(ecoCRm_lm_above1[["model"]])) {
ecoCRm_lmOutput_above1$m_slope[i] <- ecoCRm_lm_above1$model[[i]]$coefficients[2]
ecoCRm_lmOutput_above1$m_se[i] <- summary(ecoCRm_lm_above1$model[[i]])$coefficient["m", "Std. Error"]
}
ecoCRm_lmOutput_above1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,
logLik:deviance,
m_slope:m_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
m_slope = "m Coeff.",
m_se = "m SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers'
Pool (m)"),
subtitle = md("C:R values between (1,10].")) %>%
gtsave(., filename = "../Results/CR_m_lm_above1.rtf")
ecoCR_long %>%
filter(., CRratio > 1) %>%
ggplot(., aes(x = m, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values between (1,10].") +
scale_color_manual(name = "Legend", values = c("cyan", "magenta")) +
theme_pubr(legend = "right", base_size = 10)
Figure 15: Distribution of C:R (1, 10] as movement rate from the dispersers’ pool to ecosystem 2 (m) increases. Focusing on only C:R values indicative of inverted, top-heavy biomass pyramids shows weakly increasing trends in the value of C:R in ecosystem 1 when consumers move against-gradient (top-left panel), and in ecosystem 2 when consumers move along-gradient (central panel). We also note an apparent decrese in C:R values in ecosystem 1 as consumers move in an homogeneous meta-ecosystem (bottom-left panel). Note, however, that none of these trends appears to influence biomass dynamics at the meta-ecosystem scale (right-most column). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVm_binned_above1.pdf", device = "pdf",
dpi = 300, width = 12, height = 12)
And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).
ecoCRm_lm_below1 <- ecoCR_long %>%
filter(., CRratio <= 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ m, data = data)))
ecoCRm_lmOutput_below1 <- ecoCRm_lm_below1 %>% summarise(glance(model))
ecoCRm_lmOutput_below1$m_slope <- NA
ecoCRm_lmOutput_below1$m_se <- NA
for (i in 1:length(ecoCRm_lm_below1[["model"]])) {
ecoCRm_lmOutput_below1$m_slope[i] <- ecoCRm_lm_below1$model[[i]]$coefficient[2]
ecoCRm_lmOutput_below1$m_se[i] <- summary(ecoCRm_lm_below1$model[[i]])$coefficient["m", "Std. Error"]
}
ecoCRm_lmOutput_below1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,logLik:deviance, m_slope:m_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
m_slope = "m Coeff.",
m_se = "m SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers'
Pool (m)"),
subtitle = md("C:R values between [0,1].")) %>%
gtsave(., filename = "../Results/CR_m_lm_below1.rtf")
ecoCR_long %>%
filter(., CRratio <= 1) %>%
ggplot(., aes(x = m, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values between [0,1].") +
theme_pubr(legend = "right", base_size = 10)
Figure 16: Distribution of C:R [0,1] as consumer movement rate from the dispersers pool to ecosystem 2 (m) increases. Here, we notice apparently increasing trends in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). We do not see any influence of these trends at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVm_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
Here, we investigate how the \(C{:}R\) biomass ratio changes with
increasing values of c, the death rate of consumers while
in the dispersers’ pool. Consumer biomass lost while transiting through
Q is lost from the meta-ecosystem and does not contribute
to the recycling pathways at either local or regional scale. Hence, it
is important to track how biomass distribution in the system may vary
with this parameter. Here is the general relationship.
ecoCRc_lm <- ecoCR_long %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ c, data = data)))
ecoCRc_lmOutput <- ecoCRc_lm %>% summarise(glance(model))
ecoCRc_lmOutput$c_slope <- NA
ecoCRc_lmOutput$c_se <- NA
for (i in 1:length(ecoCRc_lm[["model"]])) {
ecoCRc_lmOutput$c_slope[i] <- ecoCRc_lm$model[[i]]$coefficient[2]
ecoCRc_lmOutput$c_se[i] <- summary(ecoCRc_lm$model[[i]])$coefficient["c", "Std. Error"]
}
ecoCRc_lmOutput %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,logLik:deviance, c_slope:c_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
c_slope = "c Coeff.",
c_se = "c SE") %>%
gtsave(., filename = "../Results/CR_c_lm_res.rtf")
ecoCR_long %>%
ggplot(., aes(x = c, y = CRratio, )) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Consumers death rate while in the Dispersers' Pool (c)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values limited to [0,10].") +
theme_pubr(legend = "right", base_size = 10)
Figure 17: Distribution of C:R values [0,10] as the consumer death rate in the dispersers’ pool (c) increases. No trend is apparent, aside from a weak decrease in ecosystem 2 (central and center-bottom panels) that have no influence on meta-ecosystem biomass dynamics. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVc_binned.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).
ecoCRc_lm_above1 <- ecoCR_long %>%
filter(., CRratio > 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ c, data = data)))
ecoCRc_lmOutput_above1 <- ecoCRc_lm_above1 %>% summarise(glance(model))
ecoCRc_lmOutput_above1$c_slope <- NA
ecoCRc_lmOutput_above1$c_se <- NA
for (i in 1:length(ecoCRc_lm_above1[["model"]])) {
ecoCRc_lmOutput_above1$c_slope[i] <- ecoCRc_lm_above1$model[[i]]$coefficients[2]
ecoCRc_lmOutput_above1$c_se[i] <- summary(ecoCRc_lm_above1$model[[i]])$coefficient["c", "Std. Error"]
}
ecoCRc_lmOutput_above1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,logLik:deviance, c_slope:c_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
c_slope = "c Coeff.",
c_se = "c SE") %>%
tab_header(title = md("Modeling C:R ratio vs death rate while in the
Dispersers' Pool (c)"),
subtitle = md("C:R values between (1,10].")) %>%
gtsave(., filename = "../Results/CR_c_lm_above1.rtf")
ecoCR_long %>%
filter(., CRratio > 1) %>%
ggplot(., aes(x = c, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Consumer death rate while in the Dispersers' Pool (c)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values between (1,10].") +
theme_pubr(legend = "right", base_size = 10)
Figure 18: Distribution of C:R (1,10] as consumer death rate in the dispersers’ pool (c) increases. Regardless of the type of consumer movement, we see no trend in the distribution of C:R values that identify top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVc_binned_above1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).
ecoCRc_lm_below1 <- ecoCR_long %>%
filter(., CRratio <= 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ c, data = data)))
ecoCRc_lmOutput_below1 <- ecoCRc_lm_below1 %>% summarise(glance(model))
ecoCRc_lmOutput_below1$c_slope <- NA
ecoCRc_lmOutput_below1$c_se <- NA
for (i in 1:length(ecoCRc_lm_below1[["model"]])) {
ecoCRc_lmOutput_below1$c_slope[i] <- ecoCRc_lm_below1$model[[i]]$coefficient[2]
ecoCRc_lmOutput_below1$c_se[i] <- summary(ecoCRc_lm_below1$model[[i]])$coefficient["c", "Std. Error"]
}
ecoCRc_lmOutput_below1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value,logLik:deviance, c_slope:c_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
c_slope = "c Coeff.",
c_se = "c SE") %>%
tab_header(title = md("Modeling C:R ratio vs death rate while in the
Dispersers' pool (c)"),
subtitle = md("C:R values between [0,1].")) %>%
gtsave(., filename = "../Results/CR_c_lm_below1.rtf")
ecoCR_long %>%
filter(., CRratio <= 1) %>%
ggplot(., aes(x = c, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Consumer death rate while in the Dispersers' Pool (c)") +
# ggtitle("Binned distribution of C:R",
# subtitle = "C:R values between [0,1].") +
theme_pubr(legend = "right", base_size = 10)
Figure 19: Distribution of C:R [0,1] as consumer death rate in the dispersers’ pool (c) increases. Focusing on C:R < 1 brings out a weak, negative trend in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). Neither of these trends have any influence at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
ggsave(filename = "../Results/CRratioVc_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
Here we summarize the above results and produce Figures 2, 3, and 4 shown in the main text and Figures S1 to S5 in the Supplementary Materials.
We present figures for both the raw data and the log response ratio (\(LRR_X\)) of nutrient stock and biomass accumulation, nutrient flux, and trophic compartment productivity. Our reference for calculating the \(LRR_X\) will be the results for the control scenario (Control scenario: gradient-neutral consumer movement), where environmental fertility is homogeneous between local ecosystems (i.e., I1 = I2 = 10).
Mathematically, the ratio’s formula is:
\(\displaystyle LRR_X = log_{10} \biggl(\frac{X_{i, I_1 \gtrless I_2}}{X_{i, I_1 = I_2}}\biggr)\)
where \(X \in [N, P, C]\) is the trophic compartment of interest, \(i\) is either the local or meta-ecosystem, and \(I_1 \gtrless I_2\) represents the direction of fertility gradient in the scenario of interest.
A value of \(LRR_X = 0\) indicates that the ecosystem function values in the experimental scenario does not differ from that in the control scenario. Values above 0 (white background) indicate that the ecosystem function values in the experimental scenario are larger than in the control—hence, consumer movement increases them. Viceversa, values of \(LRR_X < 0\) (grey background) indicate that ecosystem function values in the control scenario are larger than in the experimental one—hence, consumer movement decreases them.
# Here we create a custom ggplot2 theme to use in the graphs below
consmov_theme <- theme_pubr() + theme(legend.position = "bottom",
legend.text = element_text(size = 16),
legend.title = element_text(size = 17),
plot.title = element_text(size = 17),
axis.title.x = element_text(size = 17,
face = "bold"),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15),
axis.title.y = element_text(size = 17,
face = "bold"),
strip.text = element_text(size = 16))
We are working with the PRED, PREDI2, and
PREDI1 datasets we produced above (sections Control
scenario: gradient-neutral consumer movement and Experimental
scenarios: along- and against-gradient consumer movement).
First, we will check whether there is overlap among the equilibria
that are either unstable or without biological meaning across
PRED, PREDI2, and PREDI1. Then,
we will join them into a new object, EnvFert_res.
Mode FALSE TRUE
logical 61 94
Mode FALSE TRUE
logical 29 126
There is not overlap among the three datasets in terms of equilibria that are unstable and lack biological meaning. This is something that will need accounting for later on when we work with \(LRR_X\).
For now, we can remove all equilibria that are either unstable or do not have biological sense as we are going to create a figure using the raw data.
EnvFert_respos <- filter(EnvFert_res, biosense == "yes")
We now transform EnvFert_respos to a long-format
dataframe and add two categorical variables: Function and
Scale. Function tells us which ecosystem
function is measured. Scale shows instead whether the
function is being measured at the local or meta-ecosystem scale.
EnvFert_long <- select(EnvFert_respos, !u1:e2) %>%
pivot_longer(., cols = N1:PROD_Ctot,
names_to = "EcoFunction",
values_to = "Value") %>%
dplyr::mutate(., Function = if_else(EcoFunction == "N1" |
EcoFunction == "N2" |
EcoFunction == "P1" |
EcoFunction == "P2" |
EcoFunction == "C1" |
EcoFunction == "C2" |
EcoFunction == "Ntot" |
EcoFunction == "Ptot" |
EcoFunction == "Ctot",
"Stock",
if_else(EcoFunction == "FLUX_P1" |
EcoFunction == "FLUX_P2" |
EcoFunction == "FLUX_C1" |
EcoFunction == "FLUX_C2" |
EcoFunction == "FLUX_Ptot" |
EcoFunction == "FLUX_Ctot",
"Nutrient Flux",
"Trophic Productivity"))) %>%
dplyr::mutate(., Scale = if_else(EcoFunction == "N1" |
EcoFunction == "P1" |
EcoFunction == "C1" |
EcoFunction == "FLUX_P1" |
EcoFunction == "FLUX_C1" |
EcoFunction == "PROD_P1" |
EcoFunction == "PROD_C1",
"Ecosystem 1",
if_else(EcoFunction == "N2" |
EcoFunction == "P2" |
EcoFunction == "C2" |
EcoFunction == "FLUX_P2" |
EcoFunction == "FLUX_C2" |
EcoFunction == "PROD_P2" |
EcoFunction == "PROD_C2",
"Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility),
Function = factor(Function),
Scale = factor(Scale),
EcoFunction = factor(EcoFunction))
head(EnvFert_long)
# A tibble: 6 × 9
biosense Fertility TIME I1 I2 EcoFunction Value Function
<fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl> <fct>
1 yes Equal 1 10 10 N1 3.07 Stock
2 yes Equal 1 10 10 P1 1.43 Stock
3 yes Equal 1 10 10 C1 1.31 Stock
4 yes Equal 1 10 10 N2 3.26 Stock
5 yes Equal 1 10 10 P2 1.55 Stock
6 yes Equal 1 10 10 C2 2.68 Stock
# … with 1 more variable: Scale <fct>
Now that we have a dataset ready for plotting, we produce the graphs. The next code chunk produces Figure S1 in the Supplementary Materials. Note that we show the code, but not the figure itself as it is presented in the Supplementary Materials.
StockSumm <- EnvFert_long %>% filter(., Function == "Stock") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1", Nutrients = "N2", Nutrients = "Ntot", "Primary\n Producers" = "P1", "Primary\n Producers" = "P2", "Primary\n Producers" = "Ptot", Consumers = "C1", Consumers = "C2", Consumers = "Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients", "Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 10, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral", "Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Stock and Biomass") +
ylab(" ") +
xlab(" ") +
coord_cartesian(ylim = c(0, 75))
FluxSumm <- EnvFert_long %>% filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, "Primary\n Producers" = "FLUX_P1", "Primary\n Producers" = "FLUX_P2", "Primary\n Producers" = "FLUX_Ptot", Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral", "Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Flux") +
ylab(" ") +
xlab(" ") +
coord_cartesian(ylim = c(0, 350))
ProdSumm <- EnvFert_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, "Primary\n Producers" = "PROD_P1", "Primary\n Producers" = "PROD_P2", "Primary\n Producers" = "PROD_Ptot", Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral", "Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Trophic Compartment Productivity") +
ylab(" ") +
xlab("Trophic compartment") +
coord_cartesian(ylim = c(0, 500))
FuncAll <- StockSumm + FluxSumm + ProdSumm +
plot_layout(ncol = 1, nrow = 3, guides = "collect") +
plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
consmov_theme
FuncAll
ggsave(FuncAll, filename = "../Results/FunctionSummary_10kN.pdf", device = "pdf", dpi = 300, width = 11, height = 10)
Now, we produce Figures 2, 3, and 4 from the main
text. First, we are going to create new dataframes
(PREDI2rr and PREDI1rr1) that contain only the
response ratio values for nutrient stock and biomass accumulation,
nutrient flux, and trophic compartment productivity of all state
variables. These dataframes will also contain information about the
fertility conditions under which we measure each ecosystem function.
After calculating the response ratio, we will remove the unstable,
biologically not meaningful equilibria from PREDI1 and
PREDI2 before merging them in a new object,
EnvFertRR.
In calculating \(LRR_X\), the
discrepancy in the number of equilibria that are either unstable or lack
biological meaning among PRED, PREDI1, and
PREDI2 matters. To account for this, let us first create
two objects, I1NsUsrm and I2NsUsrm that
contain the unstable, nonsensical equilibria that are found in
PRED but are not included in either PREDI1 and
PREDI2. We are going to remove the ones in
PREDI1 and PREDI2 making use of the columns
biosense and stable. Conversely, we will use
I1NsUsrm and I2NsUsrm to remove instance of
the response ratio being calculated using a denominator that is from an
unstable, nonsensical equilibrium.
Then, let’s create two new objects, PREDI2rr and
PREDI1rr, to store the response ratio values calculated by
dividing ecosystem function values in PREDI2 and
PREDI1 by those in PRED, respectively. We will
then merge these two objects in a single one called
EnvFertRR, and calculated the \(C{:}R\) biomass ratio.
# need to keep raw data in the same df as the ratio to select the >0 cases
# for state variables
PREDI2rr <- PREDI2 %>%
# we are going to batch-create new response ratio versions of each variable
# we start with the "stock" group of variables
dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1, N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot, Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
# next, we work on the "flux" group of variables
dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1, FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2, FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
# finally, we work on the "productivity" group of variables
dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1, PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2, PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
# and now let's subset the dataframe to include the resp ratio columns
dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
dplyr::filter(., biosense == "yes" | stable == "yes") %>%
dplyr::filter(., !(TIME %in% I2NsUsrm$TIME))
# repeat, but for PREDI1
PREDI1rr <- PREDI1 %>%
# we are going to batch-create new response ratio versions of each variable
# we start with the "stock" group of variables
dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1, N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot, Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
# next, we work on the "flux" group of variables
dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1, FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2, FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
# finally, we work on the "productivity" group of variables
dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1, PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2, PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
# and now let's subset the dataframe to include only the resp ratio columns
dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
dplyr::filter(., biosense == "yes" | stable == "yes") %>%
dplyr::filter(., !(TIME %in% I1NsUsrm$TIME))
EnvFertRR <- bind_rows("Along" = PREDI2rr,
"Against" = PREDI1rr,
.id = "Fertility") %>%
dplyr::select(., Fertility,
biosense,
stable,
TIME:I2,
N1rr:PROD_Crr,
N1:Ctot, u1:e2) %>%
dplyr::mutate(.,
CR_Eco1 = C1/P1,
CR_Eco2 = C2/P2,
CR_MetaEco = Ctot/Ptot
)
# now build the graph, first by pivoting to longer format, then graphing
Note that object EnvFertRR contains 19301 and not 20 000
as would be expected from merging PREDI1 and
PREDI2, because of the two rounds of removal of unstable
equilibria after said merge. The first round removes those unstable,
biologically not meaningful equilibria contained in PREDI1
and PREDI2, the second uses objects I1NsUsrm
and I2NsUsrm to remove the ones that are contained in
PRED but do not overlap with those in PREDI1
or PREDI2. T
Now, we build the graphs. First, pivot the EnvFertRR
dataframe from wide to long format.
EnvFertRR_long <- select(EnvFertRR, !(N1:e2)) %>%
pivot_longer(.,
cols = N1rr:PROD_Crr,
names_to = "Compartment",
values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C1rr" |
Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" |
Compartment == "FLUX_Crr",
"Nutrient Flux",
if_else(Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" |
Compartment == "PROD_Prr" |
Compartment == "PROD_Crr",
"Trophic Productivity",
"Stock"))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" |
Compartment == "P1rr" |
Compartment == "C1rr" |
Compartment == "FLUX_P1rr" |
Compartment == "FLUX_C1rr" |
Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr",
"Ecosystem 1",
if_else(Compartment == "N2rr" |
Compartment == "P2rr" |
Compartment == "C2rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C2rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr",
"Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility),
Function = factor(Function),
Scale = factor(Scale),
Compartment = factor(Compartment)) %>%
dplyr::mutate(., biosense = fct_drop(biosense),
stable = fct_drop(stable))
In the following code chunks, we subset and then plot each ecosystem function individually against the LRR to create Figure 2, 3, and 4 in the main text. Note that we show the code here but not the figures themselves, referring the reader to the main text to view the resulting figures.
This code chunk produces Figure 2, showing changes in nutrient stocks and biomass accumulation following consumer movement.
ptHCsubset <- c("#DDAA33", "#BB5566")
ptHCsubset_fill <- c("#f8eed6", "#f1dde0")
EFstockLRR <- EnvFertRR_long %>% filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr", Nutrients = "N2rr", Nutrients = "Nrr", "Primary\n Producers" = "P1rr", "Primary\n Producers" = "P2rr", "Primary\n Producers" = "Prr", Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients", "Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Stock and Biomass") +
ylab("LRR") +
xlab("Trophic compartment") +
consmov_theme
ggsave(EFstockLRR, filename = "../Results/Stock_LRR_10kN.pdf", device = "pdf", dpi = 300, width = 10, height = 6)
The following code chunk produces Figure 3, which shows changes to nutrient flux following consumer movement.
EFfluxLRR <- EnvFertRR_long %>% filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Primary\n Producers" = "FLUX_P1rr", "Primary\n Producers" = "FLUX_P2rr", "Primary\n Producers" = "FLUX_Prr", Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Flux") +
ylab("LRR") +
xlab("Trophic compartment") +
consmov_theme
ggsave(EFfluxLRR, filename = "../Results/Flux_LRR_10kN.pdf", device = "pdf", dpi = 300, width = 10, height = 6)
And, finally, here we produce Figure 4, that shows change in the trophic compartment productivity following consumer movement.
EFprodLRR <- EnvFertRR_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Primary\n Producers" = "PROD_P1rr", "Primary\n Producers" = "PROD_P2rr", "Primary\n Producers" = "PROD_Prr", Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Trophic Compartment Productivity") +
ylab("LRR") +
xlab("Trophic compartment") +
consmov_theme
ggsave(EFprodLRR, filename = "../Results/Production_LRR_10kN.pdf", device = "pdf", dpi = 300, width = 10, height = 6)
In the two sections below, we investigate if our results are robust to changes in biomass distribution within the system caused by moving consumers and produce Figures S2 to S5 in the Supplementary Materials.
Here, we focus on those parameter sets that produce \(C{:}R < 1\). That is, on those parameter sets that produce classic biomass pyramids (McCauley et al. 2018).
As before, we begin with graphs for the untransformed data and then proceed to the \(LRR_X\) graphs.
First, using the object EnvFert_CRbelow10 created above
that exclude all \(C{:}R\) values above
10, we further exclude all \(C{:}R\)
values above 1 and create a new object,
EnvFert_CRbelow1.
Now, we lengthen this new dataset for plotting with
ggplot2.
EnvFert_CRbelow1_long <- select(EnvFert_CRbelow1, !u1:e2) %>%
pivot_longer(., cols = N1:CR_MetaEco,
names_to = "EcoFunction",
values_to = "Value") %>%
dplyr::mutate(., Function = if_else(EcoFunction == "N1" |
EcoFunction == "N2" |
EcoFunction == "P1" |
EcoFunction == "P2" |
EcoFunction == "C1" |
EcoFunction == "C2" |
EcoFunction == "Ntot" |
EcoFunction == "Ptot" |
EcoFunction == "Ctot",
"Stock",
if_else(EcoFunction == "FLUX_P1" |
EcoFunction == "FLUX_P2" |
EcoFunction == "FLUX_C1" |
EcoFunction == "FLUX_C2" |
EcoFunction == "FLUX_Ptot" |
EcoFunction == "FLUX_Ctot",
"Nutrient Flux",
ifelse(EcoFunction == "CR_Eco1" |
EcoFunction == "CR_Eco2" |
EcoFunction == "CR_MetaEco",
"C:R",
"Trophic Productivity")))) %>%
dplyr::mutate(., Scale = if_else(EcoFunction == "N1" |
EcoFunction == "P1" |
EcoFunction == "C1" |
EcoFunction == "FLUX_P1" |
EcoFunction == "FLUX_C1" |
EcoFunction == "PROD_P1" |
EcoFunction == "PROD_C1" |
EcoFunction == "CR_Eco1",
"Ecosystem 1",
if_else(EcoFunction == "N2" |
EcoFunction == "P2" |
EcoFunction == "C2" |
EcoFunction == "FLUX_P2" |
EcoFunction == "FLUX_C2" |
EcoFunction == "PROD_P2" |
EcoFunction == "PROD_C2" |
EcoFunction == "CR_Eco2",
"Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility),
Function = factor(Function),
Scale = factor(Scale),
EcoFunction = factor(EcoFunction))
head(EnvFert_CRbelow1_long)
# A tibble: 6 × 9
biosense Fertility TIME I1 I2 EcoFunction Value Function
<fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl> <fct>
1 yes Equal 2 10 10 N1 0.773 Stock
2 yes Equal 2 10 10 P1 16.7 Stock
3 yes Equal 2 10 10 C1 0.429 Stock
4 yes Equal 2 10 10 N2 0.858 Stock
5 yes Equal 2 10 10 P2 12.4 Stock
6 yes Equal 2 10 10 C2 3.20 Stock
# … with 1 more variable: Scale <fct>
And now we use the newly created dataframe
EnvFert_CRbelow1_long to produce a figure analogous to the
Figure S1 in the Supplementary Materials.
StockSummCRbelow1 <- EnvFert_CRbelow1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction,
Nutrients = "N1",
Nutrients = "N2",
Nutrients = "Ntot",
"Primary\n Producers" = "P1",
"Primary\n Producers" = "P2",
"Primary\n Producers" = "Ptot",
Consumers = "C1",
Consumers = "C2",
Consumers = "Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction,
"Nutrients",
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 10, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F,
name = "Consumer movement",
labels = c("Equal",
"Along-gradient",
"Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Stock and Biomass") +
ylab(" ") +
xlab(" ") +
coord_cartesian(ylim = c(0, 75))
FluxSummCRbelow1 <- EnvFert_CRbelow1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction,
"Primary\n Producers" = "FLUX_P1",
"Primary\n Producers" = "FLUX_P2",
"Primary\n Producers" = "FLUX_Ptot",
Consumers = "FLUX_C1",
Consumers = "FLUX_C2",
Consumers = "FLUX_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Equal",
"Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F,
name = "Consumer movement",
labels = c("Equal",
"Along-gradient",
"Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Flux") +
ylab(" ") +
xlab(" ") +
coord_cartesian(ylim = c(0, 350))
ProdSummCRbelow1 <- EnvFert_CRbelow1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction,
"Primary\n Producers" = "PROD_P1",
"Primary\n Producers" = "PROD_P2",
"Primary\n Producers" = "PROD_Ptot",
Consumers = "PROD_C1",
Consumers = "PROD_C2",
Consumers = "PROD_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Equal",
"Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F,
name = "Consumer movement",
labels = c("Equal",
"Along-gradient",
"Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Trophic Compartment Productivity") +
ylab(" ") +
xlab("Trophic Compartment") +
coord_cartesian(ylim = c(0, 500))
FuncAllCRbelow1 <- StockSummCRbelow1 + FluxSummCRbelow1 + ProdSummCRbelow1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") +
plot_annotation( tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
consmov_theme
# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering parameter sets that produce Consumer:Resource ratios below 1.",
FuncAllCRbelow1
Figure 20: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameters sets for which C:R < 1 (unstransformed data). Notably, we observe a clear quantative difference in the nutrient stock at the meta-ecosystem scale, compared with Figure S1. However, this does not appear to cause a qualitative change in the patterns of local and meta-ecosystem functioning observe, as noted in Figure S2. All specifications as in Figure S1.
ggsave(FuncAllCRbelow1, filename = "../Results/FunctionSummary_CRbelow1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)
While the quantitative changes are more noticeable here, again in the distance between the upper and lower hinges of the boxes, there does not seem to be a marked change in the qualitative pattern observed in our results for any of the ecosystem functions of interest.
Let’s proceed to work on the \(LRR_X\) values.
As we just did for the untransformed data, let’s use the previously
created object EnvFertRR to select the parameter sets
producing \(C{:}R < 1\) and store
them in a new object, EnvFertRR_CRbelow1.
EnvFertRR_CRbelow1 <- dplyr::filter(EnvFertRR, CR_Eco1 <= 1 & CR_Eco2 <= 1 & CR_MetaEco <= 1)
We now lengthen this newly created object.
EnvFertRR_CRbelow1_long <- select(EnvFertRR_CRbelow1, !(N1:e2)) %>%
pivot_longer(.,
cols = N1rr:CR_MetaEco,
names_to = "Compartment",
values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C1rr" |
Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" |
Compartment == "FLUX_Crr",
"Nutrient Flux",
if_else(Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" |
Compartment == "PROD_Prr" |
Compartment == "PROD_Crr",
"Trophic Productivity",
ifelse(Compartment == "CR_Eco1" |
Compartment == "CR_Eco2" |
Compartment == "CR_MetaEco",
"C:R",
"Stock")))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" |
Compartment == "P1rr" |
Compartment == "C1rr" |
Compartment == "FLUX_P1rr" |
Compartment == "FLUX_C1rr" |
Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr" |
Compartment == "CR_Eco1",
"Ecosystem 1",
if_else(Compartment == "N2rr" |
Compartment == "P2rr" |
Compartment == "C2rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C2rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" |
Compartment == "CR_Eco2",
"Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility),
Function = factor(Function),
Scale = factor(Scale),
Compartment = factor(Compartment)) %>%
dplyr::mutate(., biosense = fct_drop(biosense),
stable = fct_drop(stable))
And finally we produce Figure S2, which is shown in the Supplementary Materials. Accordingly, we show the code here but not the figure itself.
# we make sure the y-axis of this figure has the same range of the "all data"
# graph above by extracting the limits of the y-axis from that graph and
# using them in this one
LRRstock_ymin <- layer_scales(EFstockLRR)$y$get_limits()[1]
LRRstock_ymax <- layer_scales(EFstockLRR)$y$get_limits()[2]
EFstockLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment,
Nutrients = "N1rr",
Nutrients = "N2rr",
Nutrients = "Nrr",
"Primary\n Producers" = "P1rr",
"Primary\n Producers" = "P2rr",
"Primary\n Producers" = "Prr",
Consumers = "C1rr",
Consumers = "C2rr",
Consumers = "Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment,
"Nutrients",
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Along",
"Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0),
fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset,
name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Stock and Biomass") +
ylab("LRR") +
xlab(" ") +
coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))
# again, we grab the y-axis limits from the relevant all-data graph
LRRflux_ymin <- layer_scales(EFfluxLRR)$y$get_limits()[1]
LRRflux_ymax <- layer_scales(EFfluxLRR)$y$get_limits()[2]
EFfluxLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment,
"Primary\n Producers" = "FLUX_P1rr",
"Primary\n Producers" = "FLUX_P2rr",
"Primary\n Producers" = "FLUX_Prr",
Consumers = "FLUX_C1rr",
Consumers = "FLUX_C2rr",
Consumers = "FLUX_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Along",
"Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0),
fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset,
name = "Consumer movement",
labels = c("Along-gradient",
"Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Flux") +
ylab("LRR") +
xlab(" ") +
coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))
# again, we grab the y-axis limits from the relevant all-data graph
LRRprod_ymin <- layer_scales(EFprodLRR)$y$get_limits()[1]
LRRprod_ymax <- layer_scales(EFprodLRR)$y$get_limits()[2]
EFprodLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment,
"Primary\n Producers" = "PROD_P1rr",
"Primary\n Producers" = "PROD_P2rr",
"Primary\n Producers" = "PROD_Prr",
Consumers = "PROD_C1rr",
Consumers = "PROD_C2rr",
Consumers = "PROD_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Along",
"Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0),
fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset,
name = "Consumer movement",
labels = c("Along-gradient",
"Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Trophic Compartment Productivity") +
ylab("LRR") +
xlab("Trophic compartment") +
coord_cartesian(ylim = c(LRRprod_ymin, LRRprod_ymax))
EFallLRR_CRbelow1 <- EFstockLRR_CRbelow1 + EFfluxLRR_CRbelow1 + EFprodLRR_CRbelow1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") +
plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
consmov_theme
# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering the LRR of parameter sets that produce Consumer:Resource ratios below 1.",
# EFallLRR_CRbelow1
ggsave(EFallLRR_CRbelow1, filename = "../Results/FunctionSummary_LRR_CRbelow1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)
And here, we produce Figure S4 from the Supplementary Materials, showing median parameter values for the subset of results where \(C{:}R < 1\). Again, we show the code but not the figure itself.
# table of median parameter values
# EnvFertRR_CRbelow1 %>% select(., Fertility, u1:e2) %>% ungroup() %>% group_by(., Fertility) %>% pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>% mutate(., Parameter = as_factor(Parameter)) %>% ungroup() %>% group_by(., Fertility, Parameter) %>% summarise(., across(.cols = value, .fns = list(mean = mean, sd = sd)), .groups = "keep") %>% pivot_wider(., names_from = "Fertility", values_from = c("value_mean", "value_sd"))%>% gt(rowname_col = "Parameter", groupname_col = "Fertility") %>% fmt_number(., columns = 2:5, decimals = 2) %>% tab_spanner(label = md("Along"), columns = c("value_mean_Along", "value_sd_Along")) %>% tab_spanner(label = md("Against"), columns = c("value_mean_Against", "value_sd_Against")) %>% cols_label(value_mean_Against = "mean", value_sd_Against = "sd", value_mean_Along = "mean", value_sd_Along = "sd")
# %>% gtsave(., filename = "../Results/MeanParam_LRR_CRbelow1.rtf")
EnvFertRR_CRbelow1 %>%
select(., Fertility, u1:e2) %>%
ungroup() %>%
group_by(., Fertility) %>%
pivot_longer(., cols = -Fertility,
names_to = "Parameter",
values_to = "value") %>%
mutate(., Parameter = as_factor(Parameter)) %>%
ggplot(., aes(x = Parameter, y = value)) +
geom_boxjitter(aes(fill = Parameter,
col = Parameter),
jitter.fill = NULL, alpha = 0.1) +
scale_color_manual(values=met.brewer("Signac", 13)) +
scale_fill_manual(values=met.brewer("Signac", 13)) +
theme_pubr(legend = "none") +
ylab("Mean") +
facet_grid(.~Fertility,
labeller = as_labeller(c(`Against` = "Against-gradient Consumer movement",
`Along` = "Along-gradient Consumer movement")))
# ggsave(filename = "../Results/MeanParam_LRR_CRbelow1.pdf", dpi = 300, device = "pdf")
Here we reproduce Figures 2, 3, and 4 from the main text and Figure S1 from the Supplementary Materials using only parameter sets that produce \(C{:}R \in (1,10)\). We exclude \(C{:}R > 10\) values, as these are likely extremely rare in real-world ecosystems. We begin, as usual, with the untransformed dataframe.
First, let’s exclude from EnvFert_respos all the
parameter sets producing \(C{:}R >
10\) and those producing \(C{:}R <
1\).
EnvFert_CRabove1 <- dplyr::mutate(EnvFert_respos,
CR_Eco1 = C1/P1, # C:R for Ecosystem 1
CR_Eco2 = C2/P2, # C:R for Ecosystem 2
CR_MetaEco = Ctot/Ptot # C:R for meta-ecosystem
) %>%
dplyr::filter(., CR_Eco1 <= 10 & CR_Eco2 <= 10 & CR_MetaEco <= 10) %>%
dplyr::filter(., CR_Eco1 > 1 & CR_Eco2 > 1 & CR_MetaEco > 1)
Now we lengthen the dataset EnvFert_CRabove1 created in
the code chunk above for plotting.
EnvFert_CRabove1_long <- select(EnvFert_CRabove1, !u1:e2) %>%
pivot_longer(.,
cols = N1:CR_MetaEco,
names_to = "EcoFunction",
values_to = "Value") %>%
dplyr::mutate(., Function = if_else(EcoFunction == "N1" |
EcoFunction == "N2" |
EcoFunction == "P1" |
EcoFunction == "P2" |
EcoFunction == "C1" |
EcoFunction == "C2" |
EcoFunction == "Ntot" |
EcoFunction == "Ptot" |
EcoFunction == "Ctot",
"Stock",
if_else(EcoFunction == "FLUX_P1" |
EcoFunction == "FLUX_P2" |
EcoFunction == "FLUX_C1" |
EcoFunction == "FLUX_C2" |
EcoFunction == "FLUX_Ptot" |
EcoFunction == "FLUX_Ctot",
"Nutrient Flux",
ifelse(EcoFunction == "CR_Eco1" |
EcoFunction == "CR_Eco2" |
EcoFunction == "CR_MetaEco",
"C:R",
"Trophic Productivity")))) %>%
dplyr::mutate(., Scale = if_else(EcoFunction == "N1" |
EcoFunction == "P1" |
EcoFunction == "C1" |
EcoFunction == "FLUX_P1" |
EcoFunction == "FLUX_C1" |
EcoFunction == "PROD_P1" |
EcoFunction == "PROD_C1" |
EcoFunction == "CR_Eco1",
"Ecosystem 1",
if_else(EcoFunction == "N2" |
EcoFunction == "P2" |
EcoFunction == "C2" |
EcoFunction == "FLUX_P2" |
EcoFunction == "FLUX_C2" |
EcoFunction == "PROD_P2" |
EcoFunction == "PROD_C2" |
EcoFunction == "CR_Eco2",
"Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility),
Function = factor(Function),
Scale = factor(Scale),
EcoFunction = factor(EcoFunction))
head(EnvFert_CRabove1_long)
# A tibble: 6 × 9
biosense Fertility TIME I1 I2 EcoFunction Value Function
<fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl> <fct>
1 yes Equal 13 10 10 N1 7.29 Stock
2 yes Equal 13 10 10 P1 0.848 Stock
3 yes Equal 13 10 10 C1 3.26 Stock
4 yes Equal 13 10 10 N2 15.5 Stock
5 yes Equal 13 10 10 P2 0.894 Stock
6 yes Equal 13 10 10 C2 7.23 Stock
# … with 1 more variable: Scale <fct>
And, finally, we reproduce Figure S1 from the Supplementary Materials.
StockSummCRabove1 <- EnvFert_CRabove1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction,
Nutrients = "N1",
Nutrients = "N2",
Nutrients = "Ntot",
"Primary\n Producers" = "P1",
"Primary\n Producers" = "P2",
"Primary\n Producers" = "Ptot",
Consumers = "C1",
Consumers = "C2",
Consumers = "Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction,
"Nutrients",
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 10, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F,
name = "Consumer movement",
labels = c("Equal",
"Along-gradient",
"Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Stock and Biomass") +
ylab(" ") +
xlab(" ") +
coord_cartesian(ylim = c(0, 75))
FluxSummCRabove1 <- EnvFert_CRabove1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction,
"Primary\n Producers" = "FLUX_P1",
"Primary\n Producers" = "FLUX_P2",
"Primary\n Producers" = "FLUX_Ptot",
Consumers = "FLUX_C1",
Consumers = "FLUX_C2",
Consumers = "FLUX_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Equal",
"Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F,
name = "Consumer movement",
labels = c("Equal",
"Along-gradient",
"Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Flux") +
ylab(" ") +
xlab(" ") +
coord_cartesian(ylim = c(0, 350))
ProdSummCRabove1 <- EnvFert_CRabove1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction,
"Primary\n Producers" = "PROD_P1",
"Primary\n Producers" = "PROD_P2",
"Primary\n Producers" = "PROD_Ptot",
Consumers = "PROD_C1",
Consumers = "PROD_C2",
Consumers = "PROD_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Equal",
"Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) +
geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) +
geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F,
name = "Consumer movement",
labels = c("Equal",
"Along-gradient",
"Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Trophic Compartment Productivity") +
ylab(" ") +
xlab("Trophic compartment") +
coord_cartesian(ylim = c(0, 500))
FuncAllCRabove1 <- StockSummCRabove1 + FluxSummCRabove1 + ProdSummCRabove1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") +
plot_annotation( tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
consmov_theme
# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering parameter sets that produce Consumer:Resource ratios above 1.",
FuncAllCRabove1
Figure 21: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameter sets for which C:R > 1 and < 10 (unstransformed data). Here, biomass and nutrient stock results appear to more clearly point towards a trophic cascade at play at the meta-ecosystem scale, compared to Figure S1. Patterns of nutrient flux and trophic compartment productivity appear qualitatively unchanged. All specifications as in Figure S1.
ggsave(FuncAllCRabove1, filename = "../Results/FunctionSummary_CRabove1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)
Compared with Figure S1, we see a clear change in the pattern of biomass accumulation—consistent with the inverted, top-heavy biomass pyramid expected when \(C{:}R > 1\) (McCauley et al. 2018).
Here, we produce Figure S3 from the Supplementary Materials. We use only parameter sets that produce \(C{:}R \in (1,10)\). First, we exclude all \(C{:}R\) values below 1.
Then, we lengthen the dataset for plotting.
EnvFertRR_CRabove1_long <- select(EnvFertRR_CRabove1, !(N1:e2)) %>%
pivot_longer(.,
cols = N1rr:CR_MetaEco,
names_to = "Compartment",
values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C1rr" |
Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" |
Compartment == "FLUX_Crr",
"Nutrient Flux",
if_else(Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" |
Compartment == "PROD_Prr" |
Compartment == "PROD_Crr",
"Trophic Productivity",
ifelse(Compartment == "CR_Eco1" |
Compartment == "CR_Eco2" |
Compartment == "CR_MetaEco",
"C:R",
"Stock")))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" |
Compartment == "P1rr" |
Compartment == "C1rr" |
Compartment == "FLUX_P1rr" |
Compartment == "FLUX_C1rr" |
Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr" |
Compartment == "CR_Eco1",
"Ecosystem 1",
if_else(Compartment == "N2rr" |
Compartment == "P2rr" |
Compartment == "C2rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C2rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" |
Compartment == "CR_Eco2",
"Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility),
Function = factor(Function),
Scale = factor(Scale),
Compartment = factor(Compartment)) %>%
dplyr::mutate(., biosense = fct_drop(biosense),
stable = fct_drop(stable))
And now we produce Figure S3. Note that, as above, we do not show the figure itself here but only the code needed to produce it.
EFstockLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment,
Nutrients = "N1rr",
Nutrients = "N2rr",
Nutrients = "Nrr",
"Primary\n Producers" = "P1rr",
"Primary\n Producers" = "P2rr",
"Primary\n Producers" = "Prr",
Consumers = "C1rr",
Consumers = "C2rr",
Consumers = "Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment,
"Nutrients",
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Along",
"Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0),
fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset,
name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Stock and Biomass") +
ylab("LRR") +
xlab(" ") +
coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))
EFfluxLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment,
"Primary\n Producers" = "FLUX_P1rr",
"Primary\n Producers" = "FLUX_P2rr",
"Primary\n Producers" = "FLUX_Prr",
Consumers = "FLUX_C1rr",
Consumers = "FLUX_C2rr",
Consumers = "FLUX_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Along",
"Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0),
fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset,
name = "Consumer movement",
labels = c("Along-gradient",
"Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Nutrient Flux") +
ylab("LRR") +
xlab(" ") +
coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))
EFprodLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment,
"Primary\n Producers" = "PROD_P1rr",
"Primary\n Producers" = "PROD_P2rr",
"Primary\n Producers" = "PROD_Prr",
Consumers = "PROD_C1rr",
Consumers = "PROD_C2rr",
Consumers = "PROD_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment,
"Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility,
"Along",
"Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0),
fill = "gray95") +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
geom_boxplot(aes(color = Fertility, fill = Fertility)) +
scale_color_manual(values = ptHCsubset,
name = "Consumer movement",
labels = c("Along-gradient",
"Against-gradient")) +
scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(.~Scale, scales = "free") +
labs(title = "Trophic Compartment Productivity") +
ylab("LRR") +
xlab("Trophic compartment") +
coord_cartesian(ylim = c(LRRprod_ymin, LRRprod_ymax))
EFallLRR_CRabove1 <- EFstockLRR_CRabove1 + EFfluxLRR_CRabove1 + EFprodLRR_CRabove1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") +
plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
consmov_theme
# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering the LRR of parameter sets that produce Consumer:Resource ratios above 1 but below 10.",
# EFallLRR_CRabove1
ggsave(EFallLRR_CRabove1, filename = "../Results/FunctionSummary_LRR_CRabove1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)
As above, we now plot the median parameter values for these parameter sets, which are shown in Figure S5 in the Supplementary Materials. Accordingly, we only show the code to produce the figure and not the figure itself.
# table of median parameter sets
# EnvFertRR_CRabove1 %>% select(., Fertility, u1:e2) %>% ungroup() %>% group_by(., Fertility) %>% pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>% mutate(., Parameter = as_factor(Parameter)) %>% ungroup() %>% group_by(., Fertility, Parameter) %>% summarise(., across(.cols = value, .fns = list(mean = mean, sd = sd)), .groups = "keep") %>% pivot_wider(., names_from = "Fertility", values_from = c("value_mean", "value_sd"))%>% gt(rowname_col = "Parameter", groupname_col = "Fertility") %>% fmt_number(., columns = 2:5, decimals = 2) %>% tab_spanner(label = md("Along"), columns = c("value_mean_Along", "value_sd_Along")) %>% tab_spanner(label = md("Against"), columns = c("value_mean_Against", "value_sd_Against")) %>% cols_label(value_mean_Against = "mean", value_sd_Against = "sd", value_mean_Along = "mean", value_sd_Along = "sd")
EnvFertRR_CRabove1 %>%
select(., Fertility, u1:e2) %>%
ungroup() %>%
group_by(., Fertility) %>%
pivot_longer(., cols = -Fertility,
names_to = "Parameter",
values_to = "value") %>%
mutate(., Parameter = as_factor(Parameter)) %>%
ggplot(., aes(x = Parameter, y = value)) +
geom_boxjitter(aes(fill = Parameter, col = Parameter),
jitter.fill = NULL, alpha = 0.1) +
scale_color_manual(values=met.brewer("Signac", 13)) +
scale_fill_manual(values=met.brewer("Signac", 13)) +
theme_pubr(legend = "none") +
ylab("Mean") +
facet_grid(.~Fertility,
labeller = as_labeller(c(`Against` = "Against-gradient Consumer movement",
`Along` = "Along-gradient Consumer movement")))
# ggsave(filename = "../Results/MeanParam_LRR_CRabove1.pdf", dpi = 300, device = "pdf")
Here, we produce summary tables for the results of the model’s iterations, using the \(LRR_X\) and \(C{:}R\) values we calculated to produce the summary plots above. For \(LRR_X\), we report median values for each ecosystem function of interest in both local ecosystems and in the meta-ecosystem.
For ease of interpretation, median \(LRR_X\) values \(< 0\) will be shown in shades of magenta based on how close they are to -1 and median \(LRR_X > 0\) will be shown in shades of green based on their closeness to 1.
First, we produce dataframes to contain the median values, one for each of stock, nutrient flux, and primary and secondary productivity.
# first, compute median of log10 ratio values
EnvFertLRR_tables <- dplyr::filter(EnvFertRR_long, Ratio>=0) %>% dplyr::group_by(., Fertility, Function, Scale, Compartment) %>% dplyr::mutate(., LRR = log10(Ratio)) %>% dplyr::summarise(., min = min(LRR), median = median(LRR), max = max(LRR), .groups = "keep")
Now, we produce the tables, separating between along and against-gradient movement, as above.
EFStockTab <- EnvFertLRR_tables %>%
dplyr::ungroup() %>%
dplyr::filter(., Function == "Stock") %>%
pivot_wider(., names_from = Fertility, values_from = c(min, median, max)) %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "C1rr", "Consumers" = "C2rr", "Consumers" = "Crr", "Producers" = "P1rr", "Producers" = "P2rr", "Producers" = "Prr", "Nutrients" = "N1rr", "Nutrients" = "N2rr", "Nutrients" = "Nrr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients", "Producers", "Consumers")) %>%
gt(groupname_col = "Scale") %>%
cols_hide("Function") %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along", "median_Along", "max_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against", "median_Against", "max_Against")) %>%
cols_label(Compartment = md("**Compartment**"), "min_Along" = md("*Min.*"), "median_Along" = md("*Median*"), "max_Along" = md("*Max.*"), "min_Against" = md("*Min.*"), "median_Against" = md("*Median*"), "max_Against" = md("*Max.*")) %>%
fmt_number(., columns = c(4:9), decimals = 2) %>%
data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>% tab_header(title = md("**Change in nutrient stock and organic biomass accumulation** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."), subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>% tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
# gtsave(EFStockTab, filename = "EFbiostock.tex", path = "../Results/")
EFStockTab
| Change in nutrient stock and organic biomass accumulation in the meta-ecosystem when consumers move against or along an environmental fertility gradient. | ||||||
|---|---|---|---|---|---|---|
| Values are expressed as median LRR for 10000 iterations of the model. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
| Min. | Median | Max. | Min. | Median | Max. | |
| Ecosystem 1 | ||||||
| Consumers | −3.16 | −0.72 | −0.70 | 0.26 | 0.26 | 2.45 |
| Nutrients | −0.70 | −0.16 | 0.00 | 0.00 | 0.12 | 0.26 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | −0.80 | 0.17 | 0.92 | −1.35 | −0.27 | 0.31 |
| Nutrients | −0.68 | 0.11 | 0.26 | −0.70 | −0.14 | 0.25 |
| Producers | 0.00 | 0.04 | 1.46 | −2.72 | −0.12 | 0.00 |
| Meta-ecosystem | ||||||
| Consumers | −1.16 | 0.08 | 0.59 | −1.08 | −0.09 | 0.57 |
| Nutrients | −0.67 | 0.03 | 0.25 | −0.68 | −0.03 | 0.25 |
| Producers | 0.00 | 0.01 | 0.90 | −0.80 | −0.02 | 0.00 |
EFFluxTab <- EnvFertLRR_tables %>%
dplyr::ungroup() %>%
dplyr::filter(., Function == "Nutrient Flux") %>%
pivot_wider(., names_from = Fertility, values_from = c(min, median, max)) %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "FLUX_C1rr", "Consumers" = "FLUX_C2rr", "Consumers" = "FLUX_Crr", "Producers" = "FLUX_P1rr", "Producers" = "FLUX_P2rr", "Producers" = "FLUX_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
gt(groupname_col = "Scale") %>%
cols_hide("Function") %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along", "median_Along", "max_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against", "median_Against", "max_Against")) %>%
cols_label(Compartment = md("**Compartment**"), "min_Along" = md("*Min.*"), "median_Along" = md("*Median*"), "max_Along" = md("*Max.*"), "min_Against" = md("*Min.*"), "median_Against" = md("*Median*"), "max_Against" = md("*Max.*"))%>%
fmt_number(., columns = c(4:9), decimals = 2) %>%
data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>% tab_header(title = md("**Change in nutrient flux** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."), subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>% tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
# gtsave(EFFluxTab, filename = "EFnutflux.tex", path = "../Results/")
EFFluxTab
| Change in nutrient flux in the meta-ecosystem when consumers move against or along an environmental fertility gradient. | ||||||
|---|---|---|---|---|---|---|
| Values are expressed as median LRR for 10000 iterations of the model. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
| Min. | Median | Max. | Min. | Median | Max. | |
| Ecosystem 1 | ||||||
| Consumers | −3.16 | −0.72 | −0.70 | 0.26 | 0.26 | 2.45 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | −0.80 | 0.17 | 0.92 | −1.35 | −0.27 | 0.31 |
| Producers | 0.00 | 0.04 | 1.46 | −2.72 | −0.12 | 0.00 |
| Meta-ecosystem | ||||||
| Consumers | −0.86 | 0.07 | 0.49 | −1.06 | −0.09 | 0.57 |
| Producers | 0.00 | 0.01 | 0.90 | −0.86 | −0.02 | 0.00 |
EFProdTab <- EnvFertLRR_tables %>%
dplyr::ungroup() %>%
dplyr::filter(., Function == "Trophic Productivity") %>%
pivot_wider(., names_from = Fertility, values_from = c(min, median, max)) %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "PROD_C1rr", "Consumers" = "PROD_C2rr", "Consumers" = "PROD_Crr", "Producers" = "PROD_P1rr", "Producers" = "PROD_P2rr", "Producers" = "PROD_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
gt(groupname_col = "Scale") %>%
cols_hide("Function") %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along", "median_Along", "max_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against", "median_Against", "max_Against")) %>%
cols_label(Compartment = md("**Compartment**"), "min_Along" = md("*Min.*"), "median_Along" = md("*Median*"), "max_Along" = md("*Max.*"), "min_Against" = md("*Min.*"), "median_Against" = md("*Median*"), "max_Against" = md("*Max.*"))%>%
fmt_number(., columns = c(4:9), decimals = 2) %>%
data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>% tab_header(title = md("**Change in trophic compartment productivity** in the meta-ecosystem when consumers move _along-_ or _against-gradient_ of nutrient availability."), subtitle = md("Values are expressed as median _LRR_ for 10 000 iterations of the model.")) %>% tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
# gtsave(EFProdTab, filename = "EFprod.tex", path = "../Results/")
EFProdTab
| Change in trophic compartment productivity in the meta-ecosystem when consumers move along- or against-gradient of nutrient availability. | ||||||
|---|---|---|---|---|---|---|
| Values are expressed as median LRR for 10 000 iterations of the model. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
| Min. | Median | Max. | Min. | Median | Max. | |
| Ecosystem 1 | ||||||
| Consumers | −3.16 | −0.72 | −0.70 | 0.26 | 0.26 | 2.45 |
| Producers | −0.70 | −0.16 | 0.00 | 0.00 | 0.12 | 0.26 |
| Ecosystem 2 | ||||||
| Consumers | 0.01 | 0.22 | 1.40 | −2.63 | −0.47 | −0.01 |
| Producers | 0.00 | 0.18 | 1.40 | −2.63 | −0.34 | 0.00 |
| Meta-ecosystem | ||||||
| Consumers | −0.88 | 0.03 | 0.47 | −1.02 | −0.04 | 0.70 |
| Producers | −0.47 | 0.03 | 0.88 | −0.81 | −0.04 | 0.22 |
# first, compute median of log10 ratio values
ecoCR_g_table <- bind_rows("above 1" = ecoCRg_lmOutput_above1,
"below 1" = ecoCRg_lmOutput_below1,
.id = "CRslice") %>%
ungroup() %>%
select(., CRslice:adj.r.squared, g_slope:g_se) %>%
pivot_wider(., names_from = "CRslice",
values_from = c(r.squared, adj.r.squared, g_slope, g_se)) %>%
gt(., groupname_col = "Scale") %>%
tab_spanner(label = md("**C:R (1, 10]**"),
columns = c("r.squared_above 1",
"adj.r.squared_above 1",
"g_slope_above 1",
"g_se_above 1")) %>%
tab_spanner(label = md("**C:R [0, 1]**"),
columns = c("r.squared_below 1",
"adj.r.squared_below 1",
"g_slope_below 1",
"g_se_below 1")) %>%
fmt_number(., columns = 3:10, decimals = 3) %>%
cols_label(., "r.squared_above 1"= md("R sq."),
"adj.r.squared_above 1" = md("Adj. R sq."),
"g_slope_above 1" = md("g Coeff."),
"g_se_above 1" = md("g SE"),
"r.squared_below 1"= md("R sq."),
"adj.r.squared_below 1" = md("Adj. R sq."),
"g_slope_below 1" = md("g Coeff."),
"g_se_below 1" = md("g SE")) %>%
# data_color(., columns = c(3:10),
# colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and movement from ecosystem 1.**"), subtitle = md("The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_ values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
tab_style(style = cell_text(align = "left", size = 18),
locations = cells_title())
# gtsave(ecoCR_g_table, filename = "../Results/CR_g_summary.rtf")
ecoCR_g_table
| Relationships between biomass C:R and movement from ecosystem 1. | ||||||||
|---|---|---|---|---|---|---|---|---|
| The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0. | ||||||||
| Fertility | C:R (1, 10] | C:R [0, 1] | ||||||
| R sq. | Adj. R sq. | g Coeff. | g SE | R sq. | Adj. R sq. | g Coeff. | g SE | |
| Ecosystem 1 | ||||||||
| Against | 0.083 | 0.078 | −0.728 | 0.180 | 0.119 | 0.119 | −0.017 | 0.180 |
| Along | 0.148 | 0.147 | −0.343 | 0.022 | 0.010 | 0.010 | −0.009 | 0.022 |
| Equal | 0.177 | 0.176 | −0.537 | 0.040 | 0.039 | 0.039 | −0.016 | 0.040 |
| Ecosystem 2 | ||||||||
| Against | 0.001 | 0.001 | −0.025 | 0.016 | 0.000 | 0.000 | 0.001 | 0.016 |
| Along | 0.000 | 0.000 | 0.018 | 0.015 | 0.006 | 0.005 | 0.007 | 0.015 |
| Equal | 0.001 | 0.001 | −0.027 | 0.015 | 0.001 | 0.000 | 0.003 | 0.015 |
| Meta-ecosystem | ||||||||
| Against | 0.034 | 0.033 | −0.078 | 0.012 | 0.002 | 0.001 | −0.004 | 0.012 |
| Along | 0.048 | 0.047 | −0.100 | 0.011 | 0.002 | 0.002 | −0.004 | 0.011 |
| Equal | 0.049 | 0.048 | −0.106 | 0.012 | 0.005 | 0.005 | −0.006 | 0.012 |
# first, compute median of log10 ratio values
ecoCR_m_table <- bind_rows("above 1" = ecoCRm_lmOutput_above1,
"below 1" = ecoCRm_lmOutput_below1,
.id = "CRslice") %>%
ungroup() %>%
select(., CRslice:adj.r.squared, m_slope:m_se) %>%
pivot_wider(., names_from = "CRslice",
values_from = c(r.squared, adj.r.squared, m_slope, m_se)) %>%
gt(., groupname_col = "Scale") %>%
tab_spanner(label = md("**C:R (1, 10]**"),
columns = c("r.squared_above 1",
"adj.r.squared_above 1",
"m_slope_above 1",
"m_se_above 1")) %>%
tab_spanner(label = md("**C:R [0, 1]**"),
columns = c("r.squared_below 1",
"adj.r.squared_below 1",
"m_slope_below 1",
"m_se_below 1")) %>%
fmt_number(., columns = 3:10, decimals = 3) %>%
cols_label(., "r.squared_above 1"= md("R sq."),
"adj.r.squared_above 1" = md("Adj. R sq."),
"m_slope_above 1" = md("m Coeff."),
"m_se_above 1" = md("m SE"),
"r.squared_below 1"= md("R sq."),
"adj.r.squared_below 1" = md("Adj. R sq."),
"m_slope_below 1" = md("m Coeff."),
"m_se_below 1" = md("m SE")) %>%
# data_color(., columns = c(3:10),
# colors = col_numeric(palette = "PiYG", domain = c(-,))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and movement into ecosystem 2.**"), subtitle = md("The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_ values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
tab_style(style = cell_text(align = "left", size = 18),
locations = cells_title())
# gtsave(ecoCR_m_table, filename = "../Results/CR_m_summary.rtf")
ecoCR_m_table
| Relationships between biomass C:R and movement into ecosystem 2. | ||||||||
|---|---|---|---|---|---|---|---|---|
| The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0. | ||||||||
| Fertility | C:R (1, 10] | C:R [0, 1] | ||||||
| R sq. | Adj. R sq. | m Coeff. | m SE | R sq. | Adj. R sq. | m Coeff. | m SE | |
| Ecosystem 1 | ||||||||
| Against | 0.003 | −0.002 | 0.036 | 0.048 | 0.000 | 0.000 | 0.000 | 0.001 |
| Along | 0.000 | 0.000 | −0.014 | 0.018 | 0.001 | 0.001 | −0.002 | 0.001 |
| Equal | 0.002 | 0.001 | −0.031 | 0.026 | 0.000 | 0.000 | −0.001 | 0.001 |
| Ecosystem 2 | ||||||||
| Against | 0.000 | 0.000 | −0.008 | 0.016 | 0.000 | 0.000 | 0.000 | 0.002 |
| Along | 0.001 | 0.001 | 0.026 | 0.015 | 0.016 | 0.016 | 0.012 | 0.001 |
| Equal | 0.000 | 0.000 | −0.001 | 0.015 | 0.006 | 0.006 | 0.007 | 0.001 |
| Meta-ecosystem | ||||||||
| Against | 0.000 | −0.001 | 0.000 | 0.013 | 0.000 | 0.000 | 0.000 | 0.001 |
| Along | 0.000 | 0.000 | 0.010 | 0.011 | 0.000 | 0.000 | 0.001 | 0.001 |
| Equal | 0.001 | 0.000 | −0.014 | 0.013 | 0.000 | 0.000 | 0.001 | 0.001 |
# first, compute median of log10 ratio values
ecoCR_c_table <- bind_rows("above 1" = ecoCRc_lmOutput_above1,
"below 1" = ecoCRc_lmOutput_below1,
.id = "CRslice") %>%
ungroup() %>%
select(., CRslice:adj.r.squared, c_slope:c_se) %>%
pivot_wider(., names_from = "CRslice",
values_from = c(r.squared, adj.r.squared, c_slope, c_se)) %>%
gt(., groupname_col = "Scale") %>%
tab_spanner(label = md("**C:R (1, 10]**"),
columns = c("r.squared_above 1",
"adj.r.squared_above 1",
"c_slope_above 1",
"c_se_above 1")) %>%
tab_spanner(label = md("**C:R [0, 1]**"),
columns = c("r.squared_below 1",
"adj.r.squared_below 1",
"c_slope_below 1",
"c_se_below 1")) %>%
fmt_number(., columns = 3:10, decimals = 3) %>%
cols_label(., "r.squared_above 1"= md("R sq."),
"adj.r.squared_above 1" = md("Adj. R sq."),
"c_slope_above 1" = md("c Coeff."),
"c_se_above 1" = md("c SE"),
"r.squared_below 1"= md("R sq."),
"adj.r.squared_below 1" = md("Adj. R sq."),
"c_slope_below 1" = md("c Coeff."),
"c_se_below 1" = md("c SE")) %>%
# data_color(., columns = c(3:10),
# colors = col_numeric(palette = "PiYG", domain = c(-,))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and death rate in dispersers' pool.**"), subtitle = md("The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_ values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
tab_style(style = cell_text(align = "left", size = 18),
locations = cells_title())
gtsave(ecoCR_c_table, filename = "../Results/CR_c_summary.rtf")
ecoCR_c_table
| Relationships between biomass C:R and death rate in dispersers' pool. | ||||||||
|---|---|---|---|---|---|---|---|---|
| The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0. | ||||||||
| Fertility | C:R (1, 10] | C:R [0, 1] | ||||||
| R sq. | Adj. R sq. | c Coeff. | c SE | R sq. | Adj. R sq. | c Coeff. | c SE | |
| Ecosystem 1 | ||||||||
| Against | 0.005 | −0.001 | 0.047 | 0.050 | 0.000 | 0.000 | −0.001 | 0.001 |
| Along | 0.000 | 0.000 | −0.012 | 0.019 | 0.000 | 0.000 | 0.002 | 0.001 |
| Equal | 0.000 | −0.001 | 0.000 | 0.026 | 0.000 | 0.000 | 0.001 | 0.001 |
| Ecosystem 2 | ||||||||
| Against | 0.000 | 0.000 | 0.010 | 0.016 | 0.000 | 0.000 | 0.000 | 0.002 |
| Along | 0.000 | 0.000 | −0.004 | 0.014 | 0.002 | 0.002 | −0.004 | 0.001 |
| Equal | 0.001 | 0.001 | −0.030 | 0.015 | 0.002 | 0.002 | −0.005 | 0.002 |
| Meta-ecosystem | ||||||||
| Against | 0.003 | 0.003 | 0.025 | 0.012 | 0.000 | 0.000 | 0.000 | 0.001 |
| Along | 0.000 | −0.001 | −0.003 | 0.011 | 0.001 | 0.001 | −0.003 | 0.001 |
| Equal | 0.000 | −0.001 | 0.004 | 0.012 | 0.001 | 0.001 | −0.003 | 0.001 |
Here, we produce summary tables that show how the effects of consumer movement on local and meta-ecosystem functions vary as the biomass distribution changes in the system as a result of the consumer movement itself.
We begin by collating together the required dataframes.
Then, we lengthen the dataset as this is the preferred format for
producing tables with package gt.
CRcomparison_long <- select(CRcomparison,
!c(biosense,
stable,
CR_Eco1, # we remove the C:R values as these
CR_Eco2, # are no longer needed: the information
CR_MetaEco # is store in column "CR"
)) %>%
pivot_longer(.,
cols = N1rr:PROD_Crr,
names_to = "Compartment",
values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C1rr" |
Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" |
Compartment == "FLUX_Crr",
"Nutrient Flux",
if_else(Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" |
Compartment == "PROD_Prr" |
Compartment == "PROD_Crr",
"Trophic Productivity",
"Stock"))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" |
Compartment == "P1rr" |
Compartment == "C1rr" |
Compartment == "FLUX_P1rr" |
Compartment == "FLUX_C1rr" |
Compartment == "PROD_P1rr" |
Compartment == "PROD_C1rr" ,
"Ecosystem 1",
if_else(Compartment == "N2rr" |
Compartment == "P2rr" |
Compartment == "C2rr" |
Compartment == "FLUX_P2rr" |
Compartment == "FLUX_C2rr" |
Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" ,
"Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility),
Function = factor(Function),
Scale = factor(Scale),
Compartment = factor(Compartment))
We now remove all values of the Response Ratio that are \(< 0\) and \(log_{10}\)-transform them before calculating the median.
Finally, we produce summary tables that show the differences in median \(LRR_X\) when considering all stable parameter sets, parameter sets that produce \(C{:}R > 1\), and parameter sets that produce \(C{:}R < 1\). As before, we produce one table per ecosystem function.
CRsummStock <- CRcomparison_tables %>%
ungroup() %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "C1rr", "Consumers" = "C2rr", "Consumers" = "Crr", "Producers" = "P1rr", "Producers" = "P2rr", "Producers" = "Prr", "Nutrients" = "N1rr", "Nutrients" = "N2rr", "Nutrients" = "Nrr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients", "Producers", "Consumers")) %>%
pivot_wider(., names_from = c(CR, Fertility),
values_from = median) %>%
gt(., groupname_col = "Scale") %>%
fmt_number(., columns = 4:9, decimals = 2) %>%
tab_spanner(label = md("**Along-gradient**"),
columns = c("all_Along",
"above_1_Along",
"below_1_Along")) %>%
tab_spanner(label = md("**Against-gradient**"),
columns = c("all_Against",
"above_1_Against",
"below_1_Against")) %>%
cols_label(Compartment = md("**Compartment**"),
"all_Against" = md("*All C:R*"),
"above_1_Against" = md("*C:R > 1*"),
"below_1_Against" = md("*C:R < 1*"),
"all_Along" = md("*All C:R*"),
"above_1_Along" = md("*C:R > 1*"),
"below_1_Along" = md("*C:R < 1*")) %>%
cols_hide("Function") %>%
tab_footnote(
footnote = "Includes parameter sets producing C:R > 10.",
locations = cells_column_labels(
columns = c("all_Against", "all_Along")
)
) %>% tab_footnote(
footnote = "Does not include parameter sets producing C:R > 10.",
locations = cells_column_labels(
columns = c("above_1_Against", "above_1_Along")
)
) %>%
tab_header(title = md("**Median LRR values of biomass and nutrient stocks**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
tab_style(style = cell_text(align = "left", size = 18),
locations = cells_title())
CRsummStock
| Median LRR values of biomass and nutrient stocks | ||||||
|---|---|---|---|---|---|---|
| LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
| All C:R1 | C:R > 12 | C:R < 1 | All C:R1 | C:R > 12 | C:R < 1 | |
| Ecosystem 1 | ||||||
| Consumers | −0.72 | −0.71 | −0.72 | 0.26 | 0.26 | 0.26 |
| Nutrients | −0.16 | −0.55 | −0.15 | 0.12 | 0.20 | 0.08 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | 0.17 | 0.15 | 0.10 | −0.27 | −0.13 | −0.36 |
| Nutrients | 0.11 | 0.12 | 0.02 | −0.14 | −0.09 | −0.10 |
| Producers | 0.04 | 0.03 | 0.10 | −0.12 | −0.12 | −0.15 |
| Meta-ecosystem | ||||||
| Consumers | 0.08 | −0.15 | −0.12 | −0.09 | 0.04 | −0.09 |
| Nutrients | 0.03 | −0.18 | −0.04 | −0.03 | 0.06 | 0.00 |
| Producers | 0.01 | 0.02 | 0.05 | −0.02 | −0.05 | −0.04 |
| 1 Includes parameter sets producing C:R > 10. | ||||||
| 2 Does not include parameter sets producing C:R > 10. | ||||||
gtsave(CRsummStock, filename = "../Results/CR_summary_Stock.rtf")
CRsummFlux <- CRcomparison_tables %>%
ungroup() %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "FLUX_C1rr", "Consumers" = "FLUX_C2rr", "Consumers" = "FLUX_Crr", "Producers" = "FLUX_P1rr", "Producers" = "FLUX_P2rr", "Producers" = "FLUX_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
pivot_wider(., names_from = c(CR, Fertility),
values_from = median) %>%
gt(., groupname_col = "Scale") %>%
fmt_number(., columns = 4:9, decimals = 2) %>%
tab_spanner(label = md("**Along-gradient**"),
columns = c("all_Along",
"above_1_Along",
"below_1_Along")) %>%
tab_spanner(label = md("**Against-gradient**"),
columns = c("all_Against",
"above_1_Against",
"below_1_Against")) %>%
cols_label(Compartment = md("**Compartment**"),
"all_Against" = md("*All C:R*"),
"above_1_Against" = md("*C:R > 1*"),
"below_1_Against" = md("*C:R < 1*"),
"all_Along" = md("*All C:R*"),
"above_1_Along" = md("*C:R > 1*"),
"below_1_Along" = md("*C:R < 1*")) %>%
cols_hide("Function") %>%
tab_footnote(
footnote = "Includes parameter sets producing C:R > 10.",
locations = cells_column_labels(
columns = c("all_Against", "all_Along")
)
) %>% tab_footnote(
footnote = "Does not include parameter sets producing C:R > 10.",
locations = cells_column_labels(
columns = c("above_1_Against", "above_1_Along")
)
) %>%
tab_header(title = md("**Median LRR values of nutrient flux**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
tab_style(style = cell_text(align = "left", size = 18),
locations = cells_title())
CRsummFlux
| Median LRR values of nutrient flux | ||||||
|---|---|---|---|---|---|---|
| LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
| All C:R1 | C:R > 12 | C:R < 1 | All C:R1 | C:R > 12 | C:R < 1 | |
| Ecosystem 1 | ||||||
| Consumers | −0.72 | −0.71 | −0.72 | 0.26 | 0.26 | 0.26 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | 0.17 | 0.15 | 0.10 | −0.27 | −0.13 | −0.36 |
| Producers | 0.04 | 0.03 | 0.10 | −0.12 | −0.12 | −0.15 |
| Meta-ecosystem | ||||||
| Consumers | 0.07 | −0.01 | −0.07 | −0.09 | 0.02 | −0.11 |
| Producers | 0.01 | 0.02 | 0.04 | −0.02 | −0.05 | −0.04 |
| 1 Includes parameter sets producing C:R > 10. | ||||||
| 2 Does not include parameter sets producing C:R > 10. | ||||||
gtsave(CRsummFlux, filename = "../Results/CR_summary_Flux.rtf")
CRsummProd <- CRcomparison_tables %>%
ungroup() %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "PROD_C1rr", "Consumers" = "PROD_C2rr", "Consumers" = "PROD_Crr", "Producers" = "PROD_P1rr", "Producers" = "PROD_P2rr", "Producers" = "PROD_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
pivot_wider(., names_from = c(CR, Fertility),
values_from = median) %>%
gt(., groupname_col = "Scale") %>%
fmt_number(., columns = 4:9, decimals = 2) %>%
tab_spanner(label = md("**Along-gradient**"),
columns = c("all_Along",
"above_1_Along",
"below_1_Along")) %>%
tab_spanner(label = md("**Against-gradient**"),
columns = c("all_Against",
"above_1_Against",
"below_1_Against")) %>%
cols_label(Compartment = md("**Compartment**"),
"all_Against" = md("*All C:R*"),
"above_1_Against" = md("*C:R > 1*"),
"below_1_Against" = md("*C:R < 1*"),
"all_Along" = md("*All C:R*"),
"above_1_Along" = md("*C:R > 1*"),
"below_1_Along" = md("*C:R < 1*")) %>%
cols_hide("Function") %>%
tab_footnote(
footnote = "Includes parameter sets producing C:R > 10.",
locations = cells_column_labels(
columns = c("all_Against", "all_Along")
)
) %>% tab_footnote(
footnote = "Does not include parameter sets producing C:R > 10.",
locations = cells_column_labels(
columns = c("above_1_Against", "above_1_Along")
)
) %>%
tab_header(title = md("**Median LRR values of trophic compartment productivity**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
tab_style(style = cell_text(align = "left", size = 18),
locations = cells_title())
CRsummProd
| Median LRR values of trophic compartment productivity | ||||||
|---|---|---|---|---|---|---|
| LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
| All C:R1 | C:R > 12 | C:R < 1 | All C:R1 | C:R > 12 | C:R < 1 | |
| Ecosystem 1 | ||||||
| Consumers | −0.72 | −0.71 | −0.72 | 0.26 | 0.26 | 0.26 |
| Producers | −0.16 | −0.55 | −0.15 | 0.12 | 0.20 | 0.08 |
| Ecosystem 2 | ||||||
| Consumers | 0.22 | 0.20 | 0.22 | −0.47 | −0.28 | −0.57 |
| Producers | 0.18 | 0.17 | 0.16 | −0.34 | −0.27 | −0.34 |
| Meta-ecosystem | ||||||
| Consumers | 0.03 | −0.07 | −0.14 | −0.04 | 0.06 | −0.05 |
| Producers | 0.03 | 0.01 | 0.03 | −0.04 | −0.02 | −0.05 |
| 1 Includes parameter sets producing C:R > 10. | ||||||
| 2 Does not include parameter sets producing C:R > 10. | ||||||
gtsave(CRsummProd, filename = "../Results/CR_summary_Prod.rtf")
This notebook was compiled in 3.7 minutes.
Text and figures are licensed under Creative Commons Attribution CC BY-SA 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".